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Synopsis 
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Thesis Title: Optimum Spacing of Vertical Parallel Plate Channels 
with Natural and Mixed Convection Cooling 

Name of Thesis Supervisor; Professor V. Eswaran 

Month & Year of Submission: June, 2001 


This thesis is concerned with the optimum inter-plate spacing in vertical 
channels cooled by natural and mixed convective flows, in the particular con- 
text of applications in electronic cooling. The objective is to determine the 
optimum spacing between the channel walls that maximizes the heat trans- 
fer for various Crashof numbers (Gr^-), in natural convection, and various 
combinations of Crashof and Reynolds {Ren) numbers, in mixed convection. 

For the single channel case this is equivalent to maximising the average 
wall Nusselt number (Nuav), while in the multiple channel case it involves 
maximising the product of NuavAr, where Ar is the aspect ratio of the chan- 
nel. The computations are done on a single channel that is assumed to be one 
of an infinite array of similar channels. The single channel and multiple chan- 
nel cases are essentially the volume unconstrained and volume constrained 
optimisation of the inter-plate spacing in an array of identical channels. 
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The flow is assumed to be two dimensional, incompressible and laminar, 
with negligible viscous dissipation. The Boussinesq approximation is as- 
sumed to be valid. The approach followed is to numerically solve the stream 
function, vorticity and temperature equations in the channel for flow con- 
ditions of natural and mixed convection, at various parameter values. The 
spectral element method of Patera (1984) is applied to the governing equa- 
tions in pseudo-transient form. The latter are solved by time-stepping the 
stream function, vorticity and temperature fields to steady state, from arbi- 
trary initial conditions. The average Nusselt number is then obtained from 
the steady state temperature field, and the optimum spacing is determined 
from computations in an appropriate range of Ar. 

The optimisation for mixed convection is simpler to perform because the 
volume flow rate through the channel is known a -prion. However in the 
natural convection case, the flow rate is unknown and has to be determined 
iteratively by obtaining a solution that balances the inertial, viscous, buoy- 
ancy and pressure forces in the channel. 

The numerical simulations are done for Gth values of 10^, 10"^, and 10® 
for a Prandtl number of 0.70. For the mixed convection case, simulations are 
done for Gth = 10®, 10'^, and 10®, and buoyancy parameter (Gru/Refj) val- 
ues ranging from 0.1 to 1.5. In the natural convection multiple channel case 
and the mixed convection single channel case clearly optimal spacings are 
obtained. In the other cases either weak or asymptotic optima are observed. 
In the natural convection case, the flow rate to Ar at the various Grjj is also 
determined. The velocity and temperature fields give a graphic explanation 
of the physical processes that go into determining the opt imum spaoing. 

Chapter 1 ( “Introduction” ) of the thesis underscores the importance of the 
thesis problem in view of the great and continuing interest in electronic cool- 
ing. The latter has often been stylised, in a vast body of literature, in terms 
of the problem of convective heat transfer in parallel plate channels with and 
without protuberances. The literature on natural and mixed convection in 
parallel plate channels is surveyed in the chapter. The smaller body of work 
specifically referring to optimum inter-plate spacing is reviewed in detail. 
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These studies have been previously attempted for natural and forced convec- 
tion flows only, making a similar study of the mixed convection problem a 
logical imperative. The present thesis seeks to remedy this lacuna and also to 
investigate aspects of the natural convection problem hitherto unexamined 
in earlier studies. 

Chapter 2 (“Mathematical Formulation”) formally presents the physical 
model assumed, the idealisations made, and the resulting mathematical model 
in terms of both the primitive variable velocity and pressure as well as the 
derived stream function and vorticity formulations. The separately nondi- 
mensionalised equations for natural and mixed convection are presented. The 
computation domain to be used and the associated boundary conditions are 
discussed in detail. The inertial, viscous, buoyancy and pressure force balance 
equation that needs to be satisfied by a natural convective flow is derived. 
The optimisation criterion for the inter-plate spacing for multiple channels 
is obtained. 

Chapter 3 (“Numerical Formulation”) presents the spectral element method 
of Patera (1984) adapted to the stream function and vorticity formulation 
used in the present study. The spectral discretisation based on the Cheby- 
shev Gauss-Lobatto grid points is described. A finite-difference discretisa- 
tion of the transient term of each governing equation, along with an implicit 
representation of the Laplacian terms, and explicit representations for the 
remaining terms, is shown to yield a Helmholtz equation. This is solved by a 
Chebyshev polynomial based variational formulation that yields a system of 
linear equations for the unknown variables at the grid points of each element. 
These elemental equations are combined, and along with the boundary con- 
ditions, solved by an LU decomposition method. A grid independence study 
demonstrates the levels of accuracy reached in the simulations. 

Chapter 4 (“Results and Discussion”) presents the numerical simulations 
for the natural and mixed convection cases separately. For natural convec- 
tion, the steady state fields are obtained iteratively, for Gr^ = 10^, 10'^, and 
10®, over a range of A^, and the optimum aspect ratio (i.e. the inverse of 
the optimum spacing) is determined for each Gvh. Both constrained and 
unconstrained volume optima are determined, and are found to be strongly 



dependent on Gth- The streamlines, velocity vectors and isotherms for typi- 
cal cases are depicted and help to elucidate the physics involved. The volume 
flow rate as a function ol Gth and Ar is also determined. A similar presen- 
tation of the mixed convection results is done for Gr^j = 10^, 10^, and 10^, 
with Gth! R e\j = 0.1, 0.5, 1.0, and 1.5. The optimum spacing is found to be 
strongly related to the Gth and buoyancy parameter values. 

The conclusions of this study and the scope for future work are presented in 
Chapter 5. 
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Chapter 1 
Introduction 


1.1 B ackgr ound 

With the continuing advances in science and technology, there has been a 
corresponding increase in the number of devices involving heat transfer in 
a wide variety of industries such as aerospace, automotive, energy produc- 
tion, electronics and many others. During their normal operation, the vital 
components (e.g., nuclear fuel rods, electronic chips) of many of these devices 
(i.e., nuclear reactors, electronic equipment) dissipate energy that may result 
in the increase in their temperature beyond a certain allowable limit. For 
example, in the electronic industry, the micro-miniaturization of Very Large 
Scale Integration (VLSI) has the consequence of packaging a very large num- 
ber of components into one very small ‘chip’; and the attendant volumetric 
heat generation has risen to extremely high levels. Although the energy in- 
put in most electronic devices does not exceed a few watts, the power density 
(volumetric heat generation) is huge due to the small volume of electronic 
chips. This power density can reach megawatts per cubic meter. As a matter 
of fact, the power per unit volume that must be dissipated by modem elec- 
tronic devices is of the same order of magnitude as that of a nuclear reactor 
(cf. Refai and Yovanovich, 1992) 

Obviously, high power densities make it difficult to meet the limiting 
constraint on the maximum temperature imposed on critical components of 
such thermal systems by considerations of fatigue, melting, thermal stresses, 
and variation in material properties (Nigen and Amon, 1995). For instance. 
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the thermal design of printed circuit boards requires that the junction tem- 
perature of electronic components mounted on boards be maintained below 
125°C, while the component case temperature should be between 85 and 
100°C by reliability considerations. If the thermal control is inadequate, the 
excessive temperature will cause the electronic components and equipments 
to malfunction or bum out. Moreover, investigations, have demonstrated 
that a single electronic component operating lO^C beyond the maximum 
temperature can reduce the reliability of some systems by as much as 50% 
(cf. Lin and Hung, 1993). A single electronic chip which represents a multi- 
layer thiTi coating substrate assembly typically dissipates a heat flux of 0.1-1 
MW-rrC^ (cf. Elperin and Rudin, 1995). Since the components of this assem- 
bly are made of different materials such as ceramics, metals, semiconductors, 
the mismatch of the thermal expansion coefficients can cause high thermal 
stresses leading to the mechanical failure of the coating (Elperin and Rudin 
1995). On the other hand, the widely differing dissipation rates from differ- 
ent chips on a single circuit board results in large temperature differences 
between these chips, again causing thermal stresses. Thus, even when the 
maximum temperature attained by the hottest chip is held within acceptable 
limit both electronic and mechanical failure could occur due to large inter- 
chip temperature difference (Garimella and Schlitz, 1995). 


The problem of heat removal from heat generating thermal systems has 
been an area of extensive research and is very likely to become more im- 
portant in the future due to the ever-increasing component volumetric heat 
generation rate. It is therefore essential to innovate and develop the most 
effective and reliable cooling techniques in order to ensure long component 
life and reliable performance. On the other hand, there is another class of 
thermal systems such as solar plate collectors, compact heat exchangers etc. 
which are exclusively used for energy collection, heating and cooling fluids 
without involving heat generation. Even though the thermal design of these 
thermal systems is generally not so complex, it is important to analyse their 
thermal performance. While substantial effort have been directed towards the 
research and development of advanced heating and/or cooling technologies, 
there is still a wide scope of understanding the fluid flow and heat transfer 
characteristics of thermal systems that require high performance heat trans- 
fer using components with progr^sively lighter weights, minimum volumes 
or accommodating shapes. 



1.2. PARALLEL PLATE CHANNELS AND ELECTRONIC COOLING 3 


1.2 Parallel Plate Channels and Electronic 
Cooling 

In a large number of cases of electronic cooling the geometric configuration 
addressed is that of an array of printed circuit boards stacked either horizon- 
tally or vertically each with an assortment of mounted electronic devices. In 
the simplest case, a pair of printed circuit boards can be reasonably treated as 
a parallel plate channel formed by smooth and negligibly thin heat generating 
boards. Indeed a large number of investigations pertaining to cooling of elec- 
tronic equipments do make this simplifying approximation for the purpose 
of general paxametric studies, (Bar-Cohen and Rohsenow, 1984; Hanzawa et 
at., 1988; Anand et al., 1992; Bejan and Sciubba, 1992; Morrone et al., 1997; 
Ledezma and Bejan, 1997). 

The studies on parallel plate channels can typically be categorized in 
terms of the type of flow, thermal configuration, and the geometry involved. 
The types of flow are firee (i.e., natural) convection, mixed convection, and 
forced convection. In natural convection, the buoyancy generated in a ther- 
mal field in the presence of a gravitational field essentially induces the flow. 
Natural convection cooling is passive in that it does not require an external 
source to provide the flow. It continues to play a prominent role in thermal 
management of low power electronic equipment because of its characteris- 
tically low operating noise, low cost, ease of maintenance, simplicity and 
absence of electromagnetic interference. Moreover, it is the only mechanism 
of heat dissipation from a heat generating plate to a coolant in the event 
of power failure. Its major disadvantage is that it usually cannot provide 
highly intensive cooling. Forced convection cooling is more suited for the 
latter cases, for here the flow is driven by an external source such as a fan or 
pump, etc. Given the intensity of flow, the cooling provided can be varied, 
unlike in natural convection. However, the major disadvantage of this mode 
of heat transfer is that an external source of power is needed. Mixed con- 
vection cooling is an intermediate mechanism of heat transfer in which both 
externally imposed as well as internally induced flows operate simultaneously. 
In a vertical channel, the main flow can be either upward or downward. The 
upward forced flow is termed ‘buoyancy assisted’ flow because of the natural 
convection created by buoyancy is in the same direction as the bulk flow. In 
contrast, the downward forced flow is called ‘buoyancy opposed’ flow as its 
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direction is opposite to that of natural convection. 

The flow parameter that characterizes natural convection is the Grashof 
number, Gr. When the Grashof number is much larger than unity, the nat- 
ural convection effects are significant. In typical cases of electronic cooling 
by natural convection the Gt, based on the height of the channel, usually 
lies in the range of 10^ - 10® (see Anand et al., 1992, and Morrone et al., 
1997, for instance). In mixed convection case, the relevant flow parameter 
quantifying the ratio of natural convection to forced convection effects is the 
buoyancy parameter, Gr/Re^, where the Grashof number and the Reynolds 
number. Re, respectively represent the vigour of the natural convection and 
the forced convection effects. Physically, the magnitude of Gr/J?e^ indicates 
the relative effect of buoyancy on forced convection. The limiting value of 
Gr/Re^ -4 0 and Gr/Re^ -> cx), correspond to the forced and natural con- 
vection limit, respectively On the other hand, when Gr/Re^ is of the order 
of unity, the buoyancy effect will be comparable to the forced flow effect 
(Gebhart et al., 1988). In this thesis, we are concerned only with cases of 
free and mixed convection in vertical parallel plate channels. 

The actual geometry of printed circuit boards with protruding devices 
having irregular shapes and sizes are very rarely simulated in theoretical stud- 
ies, for not only does such geometry complicate the problem, it also makes it 
very specific to a particular situation or configuration. Thus, the vast major- 
ity of analytical, numerical and experimental studies deal with parallel plate 
channels either with smooth plane surfaces or surfaces with regularly spaced 
rectangular protrusions simulating electronic devices. Sometimes smooth 
parallel plate cha nn els are used with regularly spaced heat sources simulat- 
ing flush-mounted devices. Thus, parallel plate channels are often used to 
give a theoretical baseline in studies of electronic cooling. In this thesis, we 
are concerned only with smooth parallel plate channels. 

The thermal boundary conditions that are used to simulate electronic 
cooling in parallel plate channels are typically either uniform temperature 
or uniform heat flux; sometimes, in order to simulate flush-mounted devices, 
these boundary conditions are applied on regularly spaced intervals on the 
walls alternating with sections where an insulated (i.e., no flux) boundary 
condition is applied to simulate the gaps between the electronic components. 
Apart from this, the channels might be considered to be either (a) symmetri- 
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cally heated simulating the mounting of the electronic devices on both sides 
of the channel walls or (b) asymmetrically heated where the uniform heat 
flux (or a relatively ‘hot’ uniform temperature) is applied only to one wall 
of the channel while the other is treated as an insulated boundary (or kept 
relatively ‘cold’ at a uniform temperature). This simulates situations where 
only one side of the printed circuit board has mounted electronic devices. 
Interestingly it has been shown by Bejan and Sciubba (1992) that the type 
of boundary conditions, i.e., constant temperature or constant heat flux, has 
little effect on the optimum spacing problem, which is the focus of this the- 
sis. Thus, the present work analyses only the symmetrically heated constant 
temperature case. 

Air-cooling is recognised as an important technique in the thermal de- 
sign of electronic packages and promises to gain in importance because it is 
simple and robust, particularly with respect to safe operation (cf. Ledezma 
and Bejan, 1996). Thermal design engineers in the electronic industry are 
constantly trying to achieve even better performance with air-cooling. 


1.3 Optimum Spacing of Parallel Plate Chan- 
nels 

In the foregoing section, it has been highlighted that heat-generating compo- 
nents are typically mounted on arrays of equidistant vertical circuit boards 
resembling heated parallel plate channels. The cooling of these circuit boards 
is crucial in maintaining a safe operating temperature and peak performance. 
However, the problem confronting the thermal design engineers is to obtain 
an adequate cooling design that balances performance with manufacturing 
and operating costs. . An ill-conceived thermal design can either be ineffec- 
tive thereby threatening the life span of a system, or over-designed, adding 
unnecessary size, weight and costs. 

The fundamental problem of heat transfer augmentation, which finds an 
important application in the thermal design of electronic packages as well, 
mainly consists of maximizing the thermal conductance between the heat 
generating surface and the coolant. For the thermal design of a stack of 
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equidistant vertical heat generating circuit boards, this can be accomplished 
by suitably selecting the heat transfer mechanism, such as free, mixed or 
forced convection, and determining the optimum spacing between the circuit 
boards. The latter can be further classified as (a) volume unconstrained, 
and (b) volume constrained, optimization of inter-plate spacing. Volume un- 
constrained optimization will be applicable to optimum spacing in a single 
channel formed by a pair of parallel plates, while volume constrained optimi- 
sation applies to several plates in a fixed volume. It is important to note that 
volume constrained optimization oflters a new opportunity as in many appli- 
cations it may be possible to vary the spacing between circuit boards stacked 
in a fixed volume such that the overall thermal conductance is maximum 
while the circuit boards are maintained within safe operating temperature 
limit. 

The optimum spacing between heated parallel plates depends on a number 
of thermo-geometric parameters such as the coolant Prandtl number, height 
of the plates, thermal and geometric configuration of the plates and, most 
importantly, the Grashof number (for natural convection cooling) or both the 
Grashof number and the buoyancy parameter (for mixed convection cooling). 

In the single channel case, as the spacing between the plates is increased 
keeping other parameters fixed, the coolant volume flow rate can be expected 
to increase. This increase is more unambiguous in mixed convection cooling, 
where the flow is partially imposed, in contrast to natural convection case. 
The increase in volume flow rate manifests in greater heat transfer rate. How- 
ever, this increase in heat transfer rate is not e3q)ected to occur indefinitely 
with increase in spacing. As the spacing between the plates become large, 
the heat transfer rate from each plate tends to attain the single plate limit 
attained by a solitary plate in an infinite medium. In fact at somewhat lower 
values of inter-plate spacing the heated fluid causes upward suction popularly 
known as the ‘chimney effect’ that accelerates the flow and increase the heat 
transfer rate. Thus, the heat transfer in a parallel plate channel typically 
shows an increase and then a decrease with increasing inter-plate spacing, 
the latter decrease being not always clearly discemable. This optimum spac- 
ing is obviously of practical interest because it would give the maximum heat 
transfer in the channel. Its relevance to electronic cooling is obvious. 

In the multiple channel case, a slightly different problem needs to be ad- 
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dressed. Here we are concerned with the maximization of heat transfer from 
a series of equal-spaced parallel plates in a fixed volume. For an allowable 
maximum temperature rise, the heat transferred to the channel decreases 
with decreasing spacing, but the total number of plates increases. Hence 
there is an optimum number of plates per unit volume. Therefore, if single 
channel case is essentially an unconstrained optimisation of heat transfer in a 
channel, the multiple channel case amounts to optimisation of heat transfer 
constrained by the space being occupied. It will be shown later that in the 
single channel case, the optimum is essentially that of the maximum average 
Nusselt number, Nuav while in the multiple channel case, it is obtained at 
the maximum of the product of average Nusselt number and aspect ratio, 
i.e., NuavAr where Ar is the aspect ratio of the channel. The main thrust 
of this thesis is to determine the optimum spacing in single channel as well 
as multiple channels with natural and mixed convection cooling using air as 
the coolant. 

1.4 Previous Studies on Parallel Plate Chan- 
nels 

Convective heat transfer in heated parallel plate channels is encountered in 
a wide variety of practical applications, such as heat exchanger design, solar 
energy collection, heating of building via Trombe walls, cooling of nuclear 
reactors and electronic packages. During the past several decades, these im- 
portant apphcations — particularly, the cooling of modern electronic equip- 
ment — have motivated a vast amount of fundamental research into the 
structure of free, forced, and mixed convection flows through heated parallel 
plate channels, and of the optimum geometric arrangement of heated parallel 
plates cooled by natural and forced convection. 

Since the focus in the present thesis is on the optimum spacing between 
symmetrically heated vertical parallel plates cooled by natural and mixed 
convection, only those papers, from a vast literature on parallel plate chan- 
nels, most relevant to the present work will be presented in this section. In 
order to put the present study in perspective, the review is presented thus: 
in the first subsection, the investigations on natural convection in heated 
vertical parallel plates dealing with the fundamental aspects of the present 
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work are discussed. The next subsection focuses on the selected fundamental 
investigations on mixed convection in heated vertical parallel plates. Finally, 
in the third subsection, the complete literature dealing with the problem of 
optimum spacing between heated parallel plates is reviewed, classified fur- 
ther by the type of flow involved. 


1.4.1 Studies on Natural Convection 

Since the pioneering work of Elenbaas (1942), this topic has been studied 
extensively by many investigators. Analytical, experimental and numerical 
studies dealing with the effects of the heating mode (symmetrical or asym- 
metrical, uniform heat flux or uniform wall temperature), the inlet conditions 
and radiant heat exchange can be found in the literature. An excellent review 
of the previous work on this subject can be found in Peterson and Ortega 
(1990) and Manca, et al. (2000). 

A numerical study of transient laminar free convection between S 3 Tiimet- 
rically heated vertical parallel plates maintained at uniform temperature was 
done by Kettleborough (1972). In order to more realistically resolve the flow 
profile at the inlet of the channel, he was the first to use a computational 
domain that included, along with the actual physical domain between the 
two plates, a relatively large region upstream of the channel entrance. He 
also assumed that the velocity at the channel outlet was parallel to the walls 
and that the flow in the outflow region did not affect the flow field within the 
channel. Following a mixed stream function-vorticity and velocity formula- 
tion, the two-dimensional Navier-Stokes and energy equations incorporating 
the Boussinesq approximation were solved using the finite difference tech- 
nique. He concluded that the steady state velocity profile at the inlet of 
the channel was parabolic for low Grashof numbers, while for high Grashof 
numbers a minimum at the centre existed. His predicted results were in poor 
agreement with that of Elenbaas (1942). 

The same problem, including the entrance effect, was numerically in- 
vestigated by Nakamura et al. (1982). Using the stream function-vorticity 
formulation along with the Boussinesq approximation, the system of elliptic 
Navier-Stokes and energy equations was solved in a similar computational do- 
main as employed by Kettleborough (1972) with the difference that the wall 
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boundary condition for the stream function was obtained from the pressure 
drop along the channel walls. Special emphasis was given to the determina- 
tion of induced mass flow rate between the plates. Their results, although 
significantly different from Kettleborough (1972), were found in reasonable 
agreement with Elenbaas (1942). 

The effect of inter-plate spacing on laminar natural convection in a ver- 
tical channel bounded by an isothermal and an insulated unheated wall 
was studied both experimentally and numerically by Sparrow and Azevedo 
(1985). Consideration was given to the full range of operating conditions 
between the fully-developed limit and the single plate limit. The numeri- 
cal solution was carried out by taking account of both natural convection 
in the channel and conduction in the wall. It was found that the flat plate 
heat transfer does not form an upper bound for the channel heat transfer. 
The authors concluded that the channel heat transfer is particularly sensitive 
to changes in inter-plate spacing for narrow channels and at small wall-to- 
ambient temperature differences. 

A numerical study of transient natural convection in a symmetrically 
heated parallel plate vertical channel due to a step change in plate temper- 
ature was carried out by Chang and Lin (1989), including both the inlet 
and exit effects. Following the primitive variable approach, the full elliptic 
governing equations were solved employing a numerical scheme derived from 
the SIMPLER algorithm. The computational domain comprised the actual 
physical domain between the two plates, together with two relatively large 
regions respectively upstream and downstream of the channel. The time evo- 
lution of the flow and thermal fields and a correlation for the Nusselt number 
as a function of the Rayleigh number (Ra) and the channel aspect ratio were 
presented for air for Ra varying from 10^ to 10® and the aspect ratio from 5 to 
10. It was concluded that the inclusion of the outflow region at the channel 
exit does not have a significant effect on their predicted value of steady state 
heat transfer. 

The influence of the inlet velocity and pressure conditions on the natural 
convective flow within an asymmetrically heated channel was numerically 
studied by Chappidi and Eno (1990), for both uniform wall temperature and 
uniform heat flux heating conditions. Results were obtained by solving the 
boundary layer type governing equations in the physical domain between 



10 


CHAPTER 1. INTRODUCTION 


the plates employing a uniform as well as a parabolic velocity profile, with 
and without pressure defect, as the inlet boundary conditions. The authors 
observed that for low Grashof number flows, inlet conditions do not have a 
significant effect on the flow and thermal fields while for high Grashof number 
flows, the flow and thermal fields are sensitive to the inlet conditions. They 
concluded that the assumption of zero pressure defect at the inlet results in 
higher predictions of local and overall Nusselt number than those obtained 
by considering the pressure defect. 

Laminar free convection in channels formed between a series of vertical 
parallel plates with flush-mounted line heat sources was numerically investi- 
gated by Kim et al. (1991), for uniform heat flux conditions. The parabolized 
Navier-Stokes and energy equations were solved in an elemental channel do- 
main using an implicit finite difference scheme employing a uniform velocity 
profile at the channel inlet and ‘repeated’ boundary conditions for the walls. 
The conduction inside the plates separating the channels was taken into ac- 
count. The effect of repeated boundary conditions, and various geometric and 
thermal parameters (such as the thermal conductivity ratio of solid to air, 
the wall thickness to channel width ratio, and the modified Grashof number) 
on the mass flow rate, maximum surface temperature and average Nusselt 
number were discussed. The authors found that the repeated boundary con- 
dition reduces the maximum temperature, which occurs at the channel exit, 
and increases mass flow rate. It was also observed that as the wall conduc- 
tion increases, mass flow rate through the channel increases and maximum 
surface temperature as well as the average Nusselt number on the hot surface 
decrease. 

A numerical study of developing free convective flow between S 3 Tnmet- 
rically heated isothermal vertical parallel plates, including entrance effects, 
was carried out by Naylor et al. (1991). Incorporating the Boussinesq ap- 
proximation, the full elliptic Navier-Stokes and energy equations were solved 
using a finite element commercial code. In order to resolve the inflow bound- 
ary conditions, the channel together with a semicircular region placed up- 
stream of the channel entrance was used as computational domain. While 
the boundary conditions at the semicircular inflow boundary were obtained 
from Jeffrey- Hamel flow, fully-developed flow conditions were imposed at 
the channel exit. The comparison of their results to those by Kettleborough 
(1972) and Nakamura et al. (1982) showed a greater accordance with the 
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results of Nakamura et al. (1982). Considering the pressure defect at the 
entrance, parabolized Navier-Stokes and energy equations were also solved 
numerically using an explicit forward marching finite difference scheme. De- 
tailed comparisons of the parabolic and elliptic solutions showed that an 
elliptic solution is necessary to get accurate local quantities. However global 
quantities (such as total flow rate and average Nusselt number) predicted by 
the elliptic and parabolic solutions were in good agreement. 

With a view to enhance the understanding of flow and heat transfer in 
the low Rayleigh number regime, Martin et al. (1991) numerically studied 
the natural convection flow through a large array of equal-spaced vertical 
parallel plates, including entrance effects. Consideration was given to sym- 
metrically heated thin plates at a uniform wall temperature. Invoking the 
Boussinesq approximation, the full elliptic Navier-Stokes and energy equa- 
tions were solved employing a variant of the SIMPLE and SIMPLEC algo- 
rithms. The computational domain comprised the actual physical domain 
between the two plates forming a single channel together with an upstream 
plenum of width equal to the inter-plate spacing. The effect of the height of 
the plenum on the heat transfer and fluid flow characteristics of the array of 
channels was evaluated. It was shown that the Elenbaas(1942) asymptotic 
expression for low Rayleigh numbers would be valid as a very special limiting 
case. The authors concluded that the heat transfer at low Rayleigh numbers 
would depend on the shape and boundary conditions of the upstream plenum. 

Buoyancy-induced flow resulting from the symmetrical heating of two 
vertical parallel plates maintained at constant temperature was numerically 
investigated by Shyy et al. (1992), including entrance and exit effects. Full 
elliptic Navier-Stokes and energy equations in a nonorthogonal curvilinear 
coordinate system, with the Boussinesq approximation, were solved using 
a second-order accurate finite difference scheme with a self-adaptive grid 
technique. Computational domain consisted of the converging inflow region 
beneath the parallel surfaces, the region between them, and the diverging 
outflow region above. The results were obtained for a Prandtl number of 0.7, 
Grashof numbers ranging between 10^ to 10^ and channel aspect ratios of 1 
and 2. The authors observed that the downstream plume solution, especially 
the pressure distribution, did not match the interior flow field produced by 
the side entrainment of the wake emanating from the vertical channel, but 
this mismatch between the two solutions did not significantly affect the flow 
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field within the channel. 

The behaviour of natural convective flow in a symmetrically heated ver- 
tical, isothermal, parallel plate channel when adiabatic extensions of various 
sizes and shapes are added to the downstream part of the channel was anal- 
ysed both numerically and experimentally by Straatman et al. (1993). Pull 
elliptic Navier-Stokes and energy equations were solved using the finite ele- 
ment method. While fully developed flow conditions were used at the exit of 
extension, the inlet boundary conditions were based on Jeffrey-Hamel flow 
employed previously by Naylor et al. (1991). Experiments were conducted 
with ambient air using a Mach-Zehnder interferometer. An excellent agree- 
ment between experimental and numerical results was obtained. The results 
were presented for Pr = 0.7 over a wide range of wall heating conditions. 
The authors observed that the adiabatic extension resulted in overall heat 
transfer enhancement, in all cases studied. A single correlation accounting 
for the channel aspect ratio, expansion ratio, modified Rayleigh number and 
heated length ratio was presented. 

The effect of Prandtl number and aspect ratio on laminar natural con- 
vection in an asymmetrically heated vertical diannel was numerically inves- 
tigated by Hernandez et al. (1994), for the uniform wall temperature con- 
dition. The full elliptic Navier-Stokes and ^ergy equations incorporating 
the Boussinesq approximation, were solved using a finite-volume discretiza- 
tion technique. The computational domain comprised the region between 
the plates and a relatively large upstream region. The results were presented 
for a wide range of Prandtl numbers, aspect ratios, and modified Rayleigh 
numbers. Correlations of the average Nusselt number as a function of mod- 
ified Rayleigh number for different Prandtl numbers were proposed. It was 
found that for low values of the modified Rayleigh number, the average Nus- 
selt number does not follow the asymptotic behaviour corresponding to the 
fully developed regime unless the aspect ratio is sufficiently large. This was 
attributed to the upstream conduction effects outside the channel 

Pujii et al. (1994) numerically and experimentally investigated the lam- 
inar free convection flow in an array of asymmetrically heated vertical par- 
allel plate diannels, taking into account the conduction within the plates. 
While one side of each plate was subjected to a uniform heat flux the other 
side was cooled in the adjacent channel. Using stream function-vorticity ap- 
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proach, the full elliptic Navier-Stokes and energy equations in the fluid and 
the two-dimensional heat conduction equation in the plate were solved. A 
single heated plate along with the fluid regions in the two channels adjacent 
to the plate was taken as the computational domain. The experiments were 
carried out for an array of five vertical parallel plates. A close agreement 
between measured and computed induced flow rates was reported, except at 
aspect ratios below 8. Correlation expression for the local Nusselt number 
proposed could predict the plate temperature distributions within an error 
of ± 5%. 

Free convection heat transfer in multiple vertical channels formed by a 
series of symmetrically heated parallel plates of finite thickness and main- 
tained at uniform temperature was numerically studied by Floryan and No- 
vak (1995). Special emphasis was given to inflow regions below the inlet 
of the channels. Invoking the Boussinesq approximation, the full elliptic 
Navier-Stokes and energy equations in primitive variables were solved using 
the finite element FIDAP code. The computational domain comprised half 
of the plate, half of the flow region between the two plates, together with a 
plenum upstream of the channel having a width equal to that of the channel. 
The height of the plenum was taken as ten percent of the channel height. 
While the boundary conditions at the outflow boundary were assumed to be 
fully-developed conditions, an upward free stream flow condition with zero 
normal gradient of the vertical component of velocity was employed at the 
inflow boundary of the upstream plenum. Results were presented for aspect 
ratios ranging from 4 to 20 and Grashof numbers (based on the channel half 
width) up to 10®. The authors reported that the channels interact with each 
other and that the heat transfer can be significantly increased by placing 
channels sufficiently close to each other. They also showed that the overall 
heat transfer can be increased up to 18 percent. If only the entrance zone is 
considered the average Nusselt number may increase by up to 48 percent. 

The effect of adiabatic chimneys on convective heat transfer in a series of 
symmetrically heated vertical parallel plate channels maintained at uniform 
temperature was numerically investigated by Shahin and Floryan (1999). In- 
voking the Boussinesq approximations full elliptic Navier-Stokes and energy 
equations were solved using a finite element commercial code. The compu- 
tational domain was similar to that of Martin et al. (1991). The authors 
observed that the addition of adiabatic chimneys enhances the heat transfer 
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for all configurations considered. 

Most recently, Campo et al. (1999) nunaerically analysed the effect of 
adding insulated extensions to a symmetrically heated vertical parallel plate 
channel cooled by natural convection. The plates were subjected to uniform 
heat flux heating conditions and insulated extensions were appended either 
at the inlet or at the exit of the channel. Using the stream function-vorticity 
approach, the full elliptic Navier-Stokes and energy equations, invoking the 
Boussinesq approximation were solved by employing control volume method. 
While vorticity and energy equations were solved by implementing the alter- 
nating direction implicit (ADI) scheme with a false transient procedure, the 
stream function was obtained by the successive over relaxation method. The 
computational domain consisted of the region between the plates combined 
with two relatively large domains, one placed upstream of the channel inlet 
and the other downstream of the channel exit. The results were presented 
in terms of wall temperature profiles, induced mass flow rates, and pressure 
profiles. The authors observed that the insulated extensions appended at the 
exit of the channel resulted in a reduction of the maximum wall temperature. 
A dramatic increase in maximum wall temperature was observed when the 
insulated extensions was added to the inlet of the channel. 


1.4.2 Studies on Mixed Convection 

Laminax mixed convection in vertical parallel plate channels has received 
considerable attention owing to its numerous practical applications. A com- 
prehensive review of the literature on this topic can be found in Aung (1987), 
Incropera (1988), and Gebhart et al. (1988). In this section, the review of 
the most relevant previous studies are presented. 

Analysis of mixed convection in the entry region of a symmetrically heated 
vertical parallel plate channel was obtained analytically by Yao (1983), for 
both uniform wall temperature and uniform heat flux. The results, valid in 
the developing flow region, reveal information on the different length scales 
that distinguish various convective mechanisms. It was also shown that 
natural convection eventually becomes the dominant heat transfer mode if 
^ constant wall temperature, and if Gr^ > Re for constant heat 

flux walls. The author also conjectured that reverse flow may be present in 
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the fully developed region when the channel walls are maintained at uniform 
temperatures. 

Chow et al. (1984) numerically investigated the simultaneous effects of 
buoyancy and axial conduction on the forced convection in a short vertical 
parallel plate channel at low Peclet numbers. The plates were maintained 
at equal constant temperatures and both heating and cooling were consid- 
ered. In order to assess the effect of axial conduction in the upstream and 
downstream directions, adiabatic extensions were added to the parallel plate 
channel. Using the stream function and vorticity formulation, the full elliptic 
Navier-Stokes and energy equations were solved in the extended physical do- 
main employing a finite difference scheme. While parabolic velocity and uni- 
form temperature profiles were assumed at the inlet of the upstream insulated 
region, fully developed conditions were applied at the exit of downstream in- 
sulated section. The authors observed that the effect of axial conduction at 
low Peclet numbers were significant both in the cases of heating and cool- 
ing. It was also observed that for the case of heating, free convection effects 
counter axial conduction somewhat near the entrance to the heated section 
whereas in cooling, axial conduction is enhanced. The enhancement of heat 
transfer in heating was observed when free convection effects became signif- 
icant. 

Steady laminar developing mixed convection of air in a vertical paral- 
lel plate channel was numerically studied by Habchi and Acharya (1986). 
Consideration was given to symmetric heating with uniform temperature as 
well as asymmetric heating with one plate maintained at uniform tempera- 
ture and the other insulated. Assuming uniform velocity and temperature 
profiles at the inlet of the channel, the parabolic governing equations were 
solved using an implicit finite difference scheme. Results were presented for 
Rayleigh numbers ranging from 10^ to 10® and buoyancy parameters in the 
range of 0.1-5. It was shown that air temperature in the channel increases 
with increasing buoyancy parameter and decreasing Rayleigh number. The 
authors also showed that the local Nusselt number attains its maximum value 
near the inlet of the channel and increases with decreasing buoyancy param- 
eter. 

Developing laminar mixed convection, including flow reversal, in a heated 
vertical parallel plate channel was numerically investigated by Aung and 
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Worku (1986a). The channel walls were subjected to uniform, but not nec- 
essarily equal, temperatures. Employing uniform velocity and temperature 
profiles at the inflow boundary, the parabolic conservation equations were 
solved using a fully implicit finite difference scheme coupled viath a marching 
solution procedure. The results were obtained for a wide range of buoyancy 
parameter and wall temperature difference ratios, for Pr = 0.72. The au- 
thors showed that buoyancy forces could cause severe distortion in the veloc- 
ity profiles, especially under asymmetric heating condition. In the symmetric 
heating case the distortion disappeared and the velocity profile became fully 
developed with a parabolic shape at some distance from the channel entrance. 
Quantitative results for the Nusselt number and bulk temperature indicated 
that the buoyancy force and asymmetric heating have a profound influence 
on the heat transfer process. It was also shown that buoyancy dramatically 
increases the hydrodynamic entrance length but diminishes the thermal de- 
velopment distance. 

Aung and Worku (1986b) analytically analysed the fully developed lami- 
nar mixed convection including flow reversal in a vertical parallel plate chan- 
nel in whidi the two walls were maintained at uniform, but not necessarily 
equal, temperatures. The authors obtained an exact solution for velocity and 
temperature and presented a criterion for the occurrence of the flow reversal. 
The said criterion implied that for asymmetric heating, reverse flow occurs 
if the parameter Gr/Re exceeds a certain threshold value. 

A numerical study dealing with developing laminar mixed convection be- 
tween asymmetrically heated vertical parallel plate channels was carried out 
by Aung and Worku (1987). In this study, the channel walls were main- 
tained at uniform but unequal heat fluxes. The governing equations, inflow 
boundary conditions and numerical solution technique employed in this inves- 
tigation were similar to those used in their previous study (Aung and Worku, 
1986a) . The velocity and temperature profiles were presented for a wide range 
of buoyancy and asymmetric heating parameters. Prom the comparison be- 
tween the uniform wall temperature and uniform heat flux cases, it became 
evident that there are significant differences between the corresponding ve- 
locity and temperature profiles in mixed convection between asymmetrically 
heated plat^. It was shown that for a uniform heat flux heating condition, 
buoyancy introduces a lesser degree of skewness in the velocity distribution, 
and so flow reversal is more prone to occur in uniform wall temperature sit- 
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uation. In particular, no flow reversal was observed for Gr/Re up to 500. 

Steady laminar mixed convection of air in a vertical parallel plate channel 
including flow reversal was numerically investigated by Ingham et al. (1988a). 
The channel walls were subjected to equal constant temperature and both 
heating as well as cooling conditions were considered. Neglecting axial diffu- 
sion, the parabolic governing equations in the stream function and vorticity 
formulation were solved using a fully implicit finite difference scheme and a 
marching procedure. Emphasis was placed on flow reversal near the centre 
of the channel or adjacent to the plates. Numerical results in the form of 
velocity and temperature profiles, Nusselt numbers and friction factors were 
presented for Gv/ Re values ranging from -300 to 70. For large magnitudes 
of Gr/Re, reverse flow was observed; and heat transfer was greatly enhanced 
over the section containing reverse flow. 

Ingham et al. (1988b) numerically studied the steady laminar mixed con- 
vection in a vertical parallel plate channel subjected to unequal constant 
temperature walls. Both cases of heating and cooling were analysed. Un- 
like Ingham et al. (1988a)j emphasis was placed on flow reversals that were 
adjacent to the colder wall. Assuming the stream wise diffusion to be neg- 
ligible the solution of the boundary layer type equations were obtained by 
the numerical technique used in the previous study (Ingham et al., 1988a). 
The results showed good agreement with existing literature. A significant 
improvement in the heat transfer in the section of the channel with reverse 
flow was observed. As expected, buoyancy aiding flow enhanced the heat 
transfer, whereas buoyancy opposing flow inhibited it. 

The Resmolds number dependence of laminar mixed convective flow be- 
tween asymmetrically heated vertical parallel plates, with and without flow 
reversal, was numerically examined by Jeng et al. (1992). Uniform wall 
temperature heating conditions were considered. The full Navier-Stokes and 
energy equations in the primitive variables were solved using the SIMPLER 
algorithm on a staggered grid. The classical second-order upwind scheme 
was adopted to model the convective terms. While uniform velocity and 
temperature profiles were assumed at the inlet of the channel, first-order ex- 
trapolation was invoked at the exit of the channel. Results were presented 
for a wide range of asymmetry parameters, Reynolds numbers and the buoy- 
ancy parameter, Gr/Re. The authors found that for Re ^ 50, both flow 
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and thermal characteristics are independent of Re for fixed asymmetry and 
buoyancy parameters. 

A numerical study of laminar steady state buoyancy aided mixed convec- 
tion between a series of vertical equal-spaced parallel plates including wall 
conduction was carried out by Wastson et al. (1996). While one surface 
of each plate was subjected to a planar heat source, the opposite surface 
was exposed to the flowing fluid in the adjacent channel. Invoking Boussi- 
nesq approximations, the fuU elliptic Navier-Stokes and energy equations in 
primitive variables were solved using a finite volume technique. The outflow 
region above the channel exit was included in the computational domain. 
While uniform velocity and temperature profiles were imposed at the inflow 
boundary of the computational domain, the natural condition was used at 
the outflow boundary. Results in the form of velocity and temperature pro- 
files as well as flow structure downstream of the plate was presented for a 
range of Grashof numbers, Reynolds numbers, plate to fluid thermal con- 
ductivity and plate to channel width ratios, keeping the channel length to 
width ratio and the fluid Prandtl number fixed. It was shown that the ve- 
locity profiles within the channel are skewed substantially to the hot wall 
zs Grf Re increases and thermal conductivity ratio decreases. The interface 
temperatures both at the hot and cold surfaces of the plate were found to 
decrease with increasing Gr/Re. 


1.4.3 Studies on Optimum Spacing 

Determination of the optimum spacing between heated parallel plates is a 
fundamental problem in the thermal design of finned heat transfer surfaces 
and electronic packaging. This problem has received considerable attention in 
the past and promises to gain in importance in the foreseeable future mainly 
due to the increasing power densities in electronic components. Analytical, 
experimental and numerical studies on optimum spacing between a pair of 
heated parallel plates (single channel) and a series of heated parallel plates 
(multiple channels) dealing with different thermal configurations (symmet- 
rical or asymmetrical, uniform wall temperature or uniform wall heat flux), 
modes of cooling (natural or forced convection) and varying treatments of 
inlet and exit conditions can be found in the literature. In this section, a 
comprehensive review of the literature on optimum spacing between heated 
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parallel plates is presented, classified as natural convection and forced con- 
vection investigations. There has been no investigation reported for the cor- 
responding mixed convection problem. 


Natural Convection Investigations 

The problem of optimum spacing between heated vertical parallel plates 
cooled by natural convection has been the subject of many investigations. 
The earliest and the most widely quoted experimental study of natural con- 
vection in a vertical channel is that of Elenbaas (1942). He considered sym- 
metrically heated square plates subjected to uniform wall temperature condi- 
tions and recommended an expression for the optimum spacing for which the 
heat dissipation from the plate surfaces is maximum. Two decades later, the 
steady laminar natural convection between two vertical plates maintained at 
the same uniform temperature was first numerically investigated by Bodoia 
and Osterle (1962). Assuming uniform velocity and temperature profiles, 
and ambient pressure at the inlet of the channel, solutions of the boundary 
layer type governing equations were obtained using a finite difference scheme. 
Prom their numerical solutions these authors also derived a criterion for the 
optimum spacing based on maximum heat dissipation. Their overall heat 
transfer results and the optimum spacing criterion were in good agreement 
with the experimental observations of Elenbaas (1942). 

Using a somewhat different approach. Levy (1971) analytically derived 
a criterion for optimum spacing between heated vertical parallel plates dis- 
sipating heat by natural convection to the ambient air. Consideration was 
given to symmetrically heated plates subjected to uniform wall temperatures. 
The criterion obtained was based on minimizing the temperature difference 
between the plate and the inlet air with respect to the inter-plate spacing. 
For a given average heat transfer rate the optimisation fiinction was obtained 
by making a lumped system heat balance and using the Nusselt number cor- 
relation developed by Bodoia and Osterle (1962). He argued that for the 
condition of minimum temperature difference to occur, the boundary layers 
on the plat^ must not merge. In another study, Levy, et al (1975) exper- 
imentally verified the criterion for optimum spacing obtained by Levy (1971). 

Bar-Cohen and Rohsenow (1984) proposed an analytic optimisation tech- 
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nique for predicting the optimum plate-to-plate spacing in an array of heated 
vertical parallel plates cooled by natural convection. Consideration was given 
to plates subjected to both uniform wall temperature and uniform wall heat 
flux conditions. Also, both symmetric and asymmetric heating conditions 
were considered. The authors derived a composite expression for the convec- 
tive heat transfer coefficient by blending two asymptotes, one for the fully 
developed limit and the other for the single plate limit. An optimisation 
function was then obtained using a lumped system energy balance. Finally, 
optimum plate-to-plate spacing for each thermal configuration was obtained 
by maximizing the total heat dissipation per unit volume of the stack. 

Anand et al. (1992) have numerically studied the optimum spacing for 
free convection between a pair of heated vertical parallel plates by solving 
boundary layer type equations using an implicit finite difference scheme. 
Considering symmetrically as well as asymmetrically heated plates for both 
uniform wall temperature and uniform heat flux heating conditions, the re- 
sults were presented for asymmetric heating parameters ranging from 0 to 1 
and Grashof numbers ranging from 10 to 10®. Like Bar-Cohen and Rohsenow 
(1984), the criterion used for the optimum spacing corresponds to the max- 
imum average Nusselt number based on the channel height. However, the 
results of optimum spacing obtained by them were not in good agreement 
with those of Bar-Cohen and Rohsenow (1984). The differences between the 
two results was attributed to the differences in approach. 

Recently, Morrone et al. (1997) addressed the problem of optimum inter- 
plate spacing associated with a single heated vertical parallel plate channel 
cooled by the upflow natural convection of air. The authors considered the 
plates to be symmetrically heated with a uniform heat flux. In contrast to 
Anand et al. (1992), the authors included the effect of momentum and en- 
ergy diffusion outside of the channel by using the full elliptic conservation 
equations with an I-shaped enlarged computational domain consisting of the 
actual region between the plates and two relatively large domains placed up- 
stream of the entrance and downstream of the channel exit. Following the 
stream function-vorticity formulation, the system of conservation equations 
were solved by a control volume based finite difference method. The results 
pertaining to optimum spacing and induced mass flow rate were obtained for 
air {Pr = 0.71) and for a range of Grashof numbers between 10^ and 10®. The 
authors also proposed correlations relating the induced mass flow rate and 
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the optimum plate spacing to the thermo-geometrical parameters involved in 
the study. They concluded that, for all Gth, the optimum spacing obtained 
by Anand et al. (1992) over-predicted their results with relative differences 
ranging from 12.6% at Gr^ = 10^ to 6% at Gth = 10®. 

Quite recently, Ledezma and Bejan (1997) have presented both numerical 
and experimental results pertaining to the geometric optimisation of an as- 
sembly of staggered vertical plates installed in a given volume and cooled by 
laminar natural convection. The plates were assumed to be symmetrically 
heated and maintained at a uniform temperature. With the computational 
domain extended both at the upstream and downstream sections, the cou- 
pled governing equations were solved by the Galerkin finite element method. 
Considering the maximization of the overall thermal conductance between 
the plates and the quiescent surrounding fluid as the basis of the optimum, 
results were reported for a wide range of geometric parameters such as the 
horizontal spacing between adjacent columns of plates and the number of 
plates mounted in each vertical column. With the Prandtl number kept con- 
stant at 0.72, the Rayleigh number based on the vertical dimension of the 
assembly was varied in the range 10® to 10®. The authors showed that over- 
all thermal conductance has a distinct maximum with respect to inter-plate 
spacing and the optimum spacing decreases as the Rayleigh number increases. 
The authors also concluded that the arrangement with no stagger is the best. 


Forced Convection Investigations 

Due to increasing technological demand on heat transfer enhancement, the 
problem of obtaining the optimum spacing between heated parallel plates 
cooled by forced convection has been studied extensively in the recent past. 
The first study dealing with the problem is that of Bejan and Sciubba (1992). 
Considering a stack of smooth and negligibly thin equal-spaced heat gener- 
ating boards and using order of magnitude analysis the authors obtained 
correlations for optimum spacing between the boards such that the total 
heat transfer rate from the stack to the coolant is maximum. The imposed 
pressure difference across the stack was fixed and boards with two different 
thermal configurations were considered: with one surface either isothermal 
or with uniform heat flux, while the other surface was always taken as in- 
sulated. The accuracy of these correlations were assessed by an analysis of 
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the developing flow and heat transfer in each board-to-board channel. It was 
shown that the type of thermal boundary condition has only a minor effect 
on the optimal board-to-board spacing and the maximum heat transfer rate. 

Bejan et al. (1993) addressed the fundamental thermal design problem as- 
sociated with positioning of a heat generating board cooled by laminar forced 
convection in an insulated parallel plate channel so that the board operating 
temperature is minimum. Results were presented for both relatively good 
and bad thermal conductance of the board substrate. The optimum spacing 
between a heat generating board, maintained at either a uniform temperature 
or with uniform heat flux, and an insulated plate was reported. For a fixed 
rate of heat generation and also fixed spacing between the insulated parallel 
plates, the advantageous conditions of dividing a heat generating board into 
two or more equal-spaced boards inside the same channel were established. 

Mereu et al. (1993) developed, both analytically and numerically, a series 
of correlations for estimating the optimum spacing of heat generating boards 
mounted in a stack cooled by single phase laminar forced convection. The 
authors took the board thickness into account in the calculation of the op- 
timum board-to-board spacing and maximum heat transfer from the stack. 
Three different fiow configurations were considered for the stack: with fixed 
pressure drop, fixed mass flow rate, and fixed pumping power. Theoreti- 
cal results for optimum spacing and maximum overall thermal conductance 
between the stack and the coolant were validated by means of numerical sim- 
ulations of the complete flow and temperature fields. 

Morega and Bejan (1994) addressed the fundamental thermal design prob- 
lem of optimum spacing between parallel heat generating boards mounted 
equidistantly in a stack of specified volume and cooled by laminar forced 
convection. The optimum board-to-board spacing based on the maximiza- 
tion of the overall thermal conductance between the stack and the coolant 
was obtained by simulating the complete flow and temperature field numeri- 
cally in one of many identical dhannels. While one of the two surfaces of the 
channels was modelled as adiabatic, three different thermal configurations 
for the other surface were considered: with uniform heat flux, flush-mounted 
discrete heat sources, and protruding type discrete heat sources. Keeping the 
imposed pressure difference across the stack fixed, the results were presented 
for different board geometries and thermal configurations. Finally, the opti- 
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mum board-to-board spacing was correlated by a simple general expression 
that could apply to other packages with discrete heat sources. 

Morega et al. (1995) investigated the optimum geometric arrangement of 
a stack of parallel plates under forced free stream flow conditions by numeri- 
cally solving the flow and temperature fields. Considering uniform heat flux 
at the plates, the conservation equations of mass, momentum, and energy for 
steady, laminar, incompressible, two-dimensional flow of a Newtonian fluid 
were solved using a finite element software package. Optimum plate-to-plate 
spacing and the number of plates that could be installed in the fixed volume 
were obtained on the basis of minimizing the hot-spot temperature registered 
inside the stack. Choosing a stack with square cross-section and fixing the 
dimensionless plate thickness at 0.05, the numerical results were presented 
for Pr = 0.72 and Reynolds numbers in the range of 10^ to 10^. It was shown 
that the best way of positioning the plates relative to one another is by spac- 
ing them equidistantly, and that for a specified overall dimension of the stack, 
there is an optimal number of plates for which the hot-spot temperature is 
minimum. The authors also observed that as the Reynolds number increases, 
the hot-spot temperature and plate-to-plate spacing decrease. 

Fowler et al. (1997) have studied, both experimentally and numerically, 
the geometric optimisation of staggered parallel plates in a fixed volume with 
forced convection. Given that the maximum temperature inside the volume 
is within a certain allowable level, the geometric optimisation was based on 
the maximization of the overall thermal conductance between the given vol- 
ume and the given external flow. While the sole aim of the experimental 
study was to demonstrate the existence of an optimal spacing between two 
adjacent row of plates, the objective of the numerical study was to investi- 
gate the effect of four different geometric parameters — the spacing between 
the plates, the number of plates installed in one row, the plate swept length, 
and the degree to which the plates can be staggered — on the heat and fluid 
flow performance of the assembly. The numerical results were obtained for 
Pr = 0.72 and a wide range of Reynolds numbers using a computational do- 
main comprising an elemental channel with an upstream and a downstream 
section. These numerical results show the existence of an optimum spacing 
between two adjacent row of plates. It was observed that the optimum spac- 
ing decreases as Reynolds number increases and interestingly, that there is 
little benefit in dividing a fixed heat transfer area into smaller plates. 
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1.5 Objective of the Present Investigation 

The literature review presented in Section 1.4.3 reveals that very few in- 
vestigators to date have attempted to optimise the spacing between heated 
vertical parallel plates cooled by laminar free convection, while no study ad- 
dressing the problem for mixed convection cooling has been reported so far. 
With the exception of Anand et al. (1992), Ledezma and Bejan (1997), and 
Morrone et al.(1997), most of the other investigators obtained the optimisa- 
tion functions employing lumped system type heat balance with the natural 
heat transfer coeflScient known a priori. Anand et al. (1992) numerically 
predicted the optimum plate-to-plate spacing for single channel by solving 
the boundary layer type of conservation equations in a computational domain 
that coincided with the actual physical domain between the plates. There- 
fore, momentum and energy diffusion that occur outside of the channel were 
excluded. Moreover, the authors assumed zero pressure defect at the inlet of 
the channel, which causes an over-prediction of the average Nusselt number, 
particularly for high Grashof number flows, as observed by Chappidi and 
Eno (1990). 

Ledezma and Bejan (1997) focussed their attention on the volume con- 
strained geometric optimisation of an assembly of staggered, heated vertical 
parallel plates cooled by natural convection. Morrone et al. (1997) were the 
first to include entrance and exit effects when studying the problem of op- 
timum spacing between symmetrically heated vertical parallel plates cooled 
by laminar natural convection. However, they used only a uniform heat flux 
condition on the plates. While they solved the full elliptic Navier-Stokes 
and energy equations, and used upstream and downstream sections to more 
realistically resolve inflow and outflow boundary conditions at the inlet and 
exit of the channel, the applicability of their numerical simulations seems to 
be restricted to the single channel, and cannot be extended even to volume 
unconstrained optimisation of the inter-plate spacing of a series of heated ver- 
tical parallel plates, which is the common configuration in electronic cooling. 
Like Anand et al. (1992), the authors too have not attempted to optimise 
the spacing between a series of equal-spaced parallel plates stacked in a fixed 
volume. 
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Prom a historical point of view, the first investigation of natural convec- 
tion flow through a series of heated vertical parallel plates are those of Kim 
et al. (1991) and Martin et al. (1991). Subsequently, Fujii et al. (1994), 
Floryan and Novak (1995), and Shahin and Floryan (1999) for natural con- 
vection and Waston et al. (1996) for mixed convection, numerically analysed 
the physics of heat and fluid flow for more or less similar thermo-geometric 
configurations. Among these, only Martin et al. (1991), Floryan and Novak 
(1995), and Shahin and Floryan (1999) solved the full elliptic conservation 
equations, including entrance effects. Their computational domain comprised 
the actual physical domain between the plates, together with a plenum up- 
stream of the channel having a width equal to the inter-plate spacing. These 
investigations stressed the importance of the realistic elliptic model, and an 
extended computational domain, over the idealistic parabolic boundary-layer 
type equations. However, they left unexamined the aspects associated with 
the optimisation of the inter-plate spacing. 

The primary objective of the present investigations is to numerically study 
the optimum plate-to-plate distance between an equal-spaced series of heated 
vertical thin parallel plates subjected to natural and mixed convection modes 
of cooling. It is assumed that the plates are symmetrically heated and main- 
tained at a uniform temperature. The optimum spacing investigations for 
natural and mixed convection are presented as two separate case-studio. 
Emphasis is placed on volume unconstrained as well as volume constrained 
optimisation of inter-plate spacing. Accordingly, invoking the Boussinesq ap- 
proximation, the steady state velocity and temperature fields are obtained 
by solving the pseudo-transient form of the full elliptic conservation equa- 
tions in the stream function-vorticity formulation using the spectral element 
method as a solution strategy. The physical domain between a pair of thin 
plates combined with the upstream plenum, of width equal to the inter-plate 
spacing, is taken as. the computational domain in order to more realistically 
resolve the inflow boundary conditions. So as to determine the optimum 
plate-to-plate spacing, for each case considered extensive numerical simula- 
tions are carried out for a Prandtl number of 0.70 and a wide range of aspect 
ratios, Grashof numbers, and buoyancy parameters (for mixed convection). 

The mathematical formulations are presented in the next chapter. The 
spectral element method is elaborated upon in Chapter 3. The results are 
presented and discussed in Chapter 4. Chapter 5 summarises the major con- 
clusions of the study, and ends with a discussion of the scope of future work. 
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Chapter 2 

Mathematical Formulation 


2.1 Introduction 

Mathematical modelling is usually central to the analysis of engineering sys- 
tems, which are often very complicated. For a typical thermo-fluids system, 
this complexity arises mainly due to the time-dependent, multidimensional 
nature of the fluid flow and heat transfer and the complex supplementary 
conditions that govern these systems. In addition, the non-linearity of the 
flow equations, the variation of thermophysical properties and, sometimes, 
the contribution of radiative heat transfer makes the analysis all the more 
complicated. Consequently, a real system is often simplified to a stylised 
model resembling the original in shape, geometry and other physical charac- 
teristics in the gross features, but not in every detail. Thus, by the appli- 
cation of fundamental physical laws, and by incorporating approximations 
and idealisations, a mathematical model representing the physical situation 
is obtained. Such a mathematical model is generally amenable to numeri- 
cal simulations, which — hopefully, without involving exorbitant time and 
effort in computation. — give an adequate picture of the physics of the system. 

In this chapter, the formulation of the problem of natural and mixed con- 
vective flows through a symmetrically heated vertical parallel-plate channel 
is presented, with the view of determining the optimum spacing between the 
plates. The approximations and idealisations made, and their validity, are 
discussed. The mathematical model in terms of both the primitive variables 
as w^ell as stream-function and vorticity formulations are presented below for 
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the sake of completeness, although only the latter formulation is used in the 
numerical simulations. The use of the simulations results in determining the 
optimum spacing is discussed. 


2.2 Physical Model 

2.2.1 The Problem of Many Channels 

Recall from the first chapter that the problem being analysed in this thesis 
concerns the optimum spacing of vertical plates (say, printed circuit boards 
[RGBs] in electronic equipment) in order to maximise heat-transfer. In other 
words, if we had a set of RGBs how would we space them so as to maximise 
the over-all heat transfer rate for a given volume? 

To study such a general problem a number of idealisation have to be 
made. The immediate ones we assume are that 

1. the vertical plates are identical 

2. the (electronic) components actually dissipating the heat are small com- 
pared to the dimensions of the channel 

3. the transverse dimension of the plates are large enough that the flow 
may be treated as two-dimensional(2-D). 

Thus the problem actually analysed in this thesis is that of 2-D flow in 
a set of identical and parallel channels, in which the plates axe treated as 
flat, and the electronic heating is replaced by a condition of uniform wall 
temperature. 

As the channels are identical, there is no need to study more than one 
channel, for the effect of the others is obtained by the use of symmetry 
boundary conditions. Such an approach, of course, implicitly assumes that 
the number of plates are large, indeed infinite. But this assumption is justi- 
fied by the very great difficulties that would be introduced were we to assume 
otherwise, and by the loss of generality that would also follow were we to as- 
sume any other specific case. 
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Thus, now we present the problem for a single channel, which is the basic 
component of our analysis. Knowing the heat transfer rate in this single 
channel, we merely multiply by the number of channels in a given space to 
obtain the total heat transfer rate in that space. 


2.2.2 The Problem in a Single Channel 

The geometric configuration and the co-ordinate system of the vertical parallel- 
plate channel is schematically depicted in Figure 2.1. The height of the 
channel is H and the spacing between the plates forming the channel is L. 
The thickness of the plates is taken as negligible compared to the inter-plate 
spacing. The origin of the co-ordinate system is affixed at the top of the left 
plate, with x as the downward vertical co-ordinate and y as the horizontal 
co-ordinate extending transversely across the channel. The force of gravity 
thus acts along the positive x-direction. The dimension of the channel normal 
to the x — y plane is taken as large in comparison with its other dimensions, 
so the fiow is treated as two-dimensional. Both the plates of the channel 
are symmetrically heated and are maintained at a constant temperature, Tu, , 
that is higher than the ambient temperature. Too. The bottom and the top 
of the channel are exposed to the ambient air. 

The temperature gradient resulting from the heating of air within the 
channel causes a density gradient, which in turn gives rise to buoyancy forces. 
In the absence of an externally imposed pressure difference, the flow is in- 
duced entirely by this buoyancy force. As a result, the heated air in the 
channel is continuously driven upward by buoyancy and is finally discharged 
into the ambient air through the upper end of the channel. Meanwhile, the 
air from below is drawn into the channel through the lower end by the pres- 
sure defect created by the lower density of hot air within the channel. When 
the flow is entirely due to the buoyancy effect, the heat transfer is due to 
natural or free convection. However, in the pr^ence of a fan blowing air 
upwards into the channel, the flow in the channel may be both due to the 
forced flow and the buoyancy effect. In this case, it is called mixed convec- 
tion. When the buoyancy effect is negligible compared to the forced flow, 
the flow condition is called forced convection. 
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Healed air 



Figure 2.1: The Physical Model and the Coordinate System 
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2.3 Approximations and Idealisations 

The physical model described in the preceding section is a simplified model, 
with respect to the surface geometry and heating condition of the parallel 
plates, when compared with a typical channel with electronic components 
actually encountered in applications. The further approximations and ideal- 
isations made for the present investigations are that: 


1. The fluid is Newtonian, viscous, with constant thermo-physical prop- 
erties. 

2. The flow is steady, laminar and two-dimensional. 

3. The Boussinesq approximation is valid. 

4. Viscous dissipation and compression work are negligibly small. 

5. Radiative heat transfer is insignificant. 


2.3.1 Boussinesq Approximation 

In most of the analytical as well as numerical studies of natural/mixed con- 
vective flows, the classical Boussinesq approximation is customarily invoked. 
This approximation has the advantage that it bypasses the calculation of the 
density field. While all other thermophysical properties of the fluid are as- 
sumed to remain constant, the density, associated with the buoyancy force, 
is assumed to vary linearly with temperature as per the following relation: 

p„ = p[l + /3(T-r„)] (2.1) 

where ^ is the coefficient of volumetric expansion at constant pressure, Too 
is the temperature of the ambient fluid and Poo is the corresponding density. 


The use of the Boussinesq approximation is only valid when the over- 
all temperature difference within the physical domain is sufficiently small. 
Zhong et al. (1985) performed a quantitative assessment of this requirement 
for natural convection in a differentially heated square cavity. They con- 


cluded that when the non-dimensional temperature difference, 5 
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is less than 0.1, the Boussinesq approximation yields an accurate prediction 
of the flow and thermal fields. They also suggested that reasonably accurate 
values of average heat transfer rates could be obtained by using the Boussi- 
nesq approximation for <5 values up to 0.2. 

With these limits of validity, the Boussinesq approximation is probably 
valid for many applications of electronic cooling. Typically, in most modern 
electronic equipment, the maximum temperature rise should not exceed 50° 
K. Therefore, for flow through channels of circuit boards, the overall temper- 
ature difference, for an ambient temperature Too — 300 K, S wall be within 
0.2. Most investigators in their studies of natural convection between vertical 
parallel-plates have invoked the Boussinesq approximation. 


2.4 Mathematical Model 

A mathematical model comprises quantitative statements of the relationships 
among the dependent and independent variables and the relevant parame- 
ters that describe some physical phenomena (Middleman, 1998). Typically, 
a mathematical model consists of equations that govern the behaviour of the 
physical system, and the associated boundary conditions. 

Employing the approximations and the idealisations listed in Section 2.3, 
the physical model described in Section 2.2 is simulated by an equivalent 
mathematical model involving the conservation of mass, momentum, and 
energy, with appropriate boundary conditions. The mathematical model 
comprising these partial differential equations, along with their boundary 
conditions, is presented in the following sections. 


2.4.1 Governing Equations 

The governing equations describing buoyancy-driven flows and transport 
have been reviewed by Jaluria (1980), Gebhart et al. (1988) and Bejan 
(1995). For the steady, laminar, two-dimensional natural/mixed convective 
flows of a Newtonian fluid, using the Boussinesq approximation, the gov- 
erning equations of mass, momentum and energy conservation are (Ledezma 
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and Bejan, 1997): 


mass: 


rc-momentum: 


u 


du 

dx 



y-momentum: 


du dv 
dx dy 


Idp / d'^u 

pdx'^^ (9y2 / 


-g.6{T-T^) 


dv dv Idp { d'^v d^v\ 

energy: 

dT dT f^T d‘^T\ 

^ dx dy \ dx^ dy^ ) 


(2.2) 


(2.3) 


(2.4) 


(2.5) 


where (u, v) are the velocity components in the x and y directions, and p and 
T are the density and temperature, while v and a are the kinematic viscosity 
and thermal diffusivity of the fluid. The quantity p is the modified pressure 
(= Pthermodynamic— Poogx) which accounts for the hydrostatic pressure, so that 
its value at the far fleld is always constant (poo)? independent of x. It is to 
be noted that the gravitational acceleration, g, acts in the x direction, in our 
chosen convention. 


The model equations (2.2) to (2.5) form a set of coupled partial differen- 
tial equations that are elliptic in nature and require that boundary conditions 
to be imposed all around the domain boundary. 


2.4.2 Boundary Conditions 

The mathematical formulation of the problem is incomplete without pre- 
scribing a speciflc number of appropriate boundary conditions that allow 
only unique solutions of the governing equations to exist. The specification 
of mathematically correct boundary conditions that ensure the uniqueness of 
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the solution, while being compatible with the physics at the boundaries, is 
not always straightforward. In this section, the computational domain and 
the corresponding boundary conditions to be specified for the primitive vari- 
ables formulation of the physical model described in Section 2.2 are discussed. 

Before arriving at the boundary conditions at various boundaries, we 
have to first identify the solution/computational domain of the problem. 
The physical domain and the corresponding computational domain of the 
problem usually differ in one respect or the other. While the computational 
domain largely depends on the geometry of the physical domain, it may be 
extended beyond it to facilitate the specification of physically and mathe- 
matically correct boundary conditions. 

Figure 2.2 depicts the computational domain pertaining to natural/mixed 
convective flows through a symmetrically heated parallel plate channel show- 
ing all the boundaries of interest. It consists of two impermeable and ‘no- 
slip’ solid boundaries represented by the lines AB and EF and four ‘free’, i.e., 
open, boundaries represented by dotted lines BC, CD, DE and AF. While 
the open boundaries CD and AF pertain to inflow and outflow respectively, 
the open boundaries BC and DE correspond to planes of symmetry which 
exist in the plenum upstream of the channel due to presence of neighbouring 
plates, as in the case of a series of parallel plate channels. Recall that the 
solution domain is symmetrical about these axes with regards to both heat 
and fluid flow. Therefore, the axes of symmetry may be used as boundaries. 

Note that the computational domain has been extended downwards to 
include the plenum upstream of the channel. This has been done mainly to 
impose realistic inflow conditions. The extended domain has the same width 
as the channel — which appears to be appropriate for natural/ mixed convec- 
tion flows through a series of parallel plate channels. A similar approach was 
adopted by Martin et al. (1991) and Floryan and Novak (1995), and more 
recently by Shahin and Floryan (1999) in their study of natural convection 
through vertical plates. In contrast to this, Kettleborough (1972), Naylor et 
al. (1991) and, more recently, Morrone et al. (1997) have taken an enlarged 
domain near the inlet of the channel which seems to be more appropriate for 
prescribing the inflow boundary conditions for a single parallel plate channel. 
The a-ppropriate height of the extended region was determined by numerical 
experimentation, in this study, to be 10% of the channel height. Extending 
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Figure 2.2: The Computational Domain and Boundary Conditions 
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the region beyond this figure does not have any further affect on the flow 
physics, for the chosen boundary conditions. The downstream region be- 
yond the exit end of the channel has been excluded from the solution domain 
mainly due to computational economy and simplicity. 

The velocity and temperature boundary conditions at the inflow and out- 
flow boundaries are not obvious. The literature reveals that researchers have 
resorted to various idealised conditions on these boundaries. As seen below, 
the boundary conditions on all except the inflow boundary are the same for 
both natural and mixed convective flows. The inflow boundary conditions, 
however, differ with the type of flow. 

Since the solid boundaries AB and EF represent the impermeable, no-slip 
and isothermal walls of the parallel plate channel, the components of velocity 
are set to zero and the temperature is fixed at T^ at these boundaries. The 
normal gradient of pressure at the solid boundary is assumed to be zero, in 
keeping with boundary layer theory. 

The boundaries BC and DE correspond to lines of symmetry. As the 
flow across a symmetry line is zero, the component of velocity normal to the 
boundary is set to zero. Also, due to symmetry, the shear stress and the 
normal gradient of temperature and pressure on these boundaries are taken 
to be zero. At the outflow boundary, AF, the two components of velocity and 
the temperature are not known. Therefore, the flow is assumed to be fully- 
developed here and the normal gradients of the components of velocity and 
the temperature are taken to be zero. The pressure at the outflow boundary 
is taken equal to that of ambient fluid. Although these outflow conditions, 
which could affect the computed flow in the channel, are somewhat adhoc, 
they are similar to the idealisations made by previous researchers. 

The conditions on the inflow boundary CD are specified differently de- 
pending on the type of flow. For natural convection flows, the v component 
of velocity and the normal gradient of the u component of velocity are set 
to zero. The flow is assumed to enter the inflow boundary at ambient tem- 
perature and pressure. For mixed convective flows, the fluid is assumed to 
cross the inflow boundary at ambient temperature and with u=0, but the u 
component of velocity is taken to be u = —Uq (recall u is positive downwards, 
in our convention). The normal gradient of pressure is assumed to be zero. 
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An Alternate Computational Domain 

Figure 2.3 illustrates an alternate computational domain along with the ap- 
propriate boundary conditions. Note that the computational domain shown 
in Figure 2.3 is exactly a symmetric one-half of the solution domain shown 
in Figure 2.2. Here the fact that the flow domain is perfectly symmetrical 
has been used to reduce the computational domain — and the computations 
needed for solutions — by one-half. By S 3 Tnmetry, the y-derivatives of u, v, 
and T, and the y-component of velocity (u) at the physical centre-line, shown 
here as the GH boundary, will be zero. This fact is used to prescribe the 
boundary conditions shown in the figure. 

However, a danger inherent in the use of ssunmetry conditions is that 
symmetry may be violated by the flow, even if the domain is symmetric. 
This generally occurs when the flow becomes unsteady, as happens at high 
Grashof numbers. The unsteady flow-field generally is not symmetric. Too 
many symmetry boundary conditions and the use of a steady-state solver 
may thus result in pseudo steady-state solutions in cases when the flow should 
have physically become unsteady. This possibility can be guarded against 
by (a) doing simulations on the full physical domain, and (b) using a pseudo 
transient algorithm that would alert us about the onset of unsteadiness. Thus 
the computations in this thesis were done with a pseudo-transient method 
(see below), and the computational domain used was the full one shown in 
Figure 2.2 rather than the truncated one in Figure 2.3. 


2.5 Dimensionless Representation 

The mathematical model is generally cast in dimensionless form, and the 
flow variables are normalised so that their values are of the order of unity. 
The nondimensionlisation reduces the number of independent parameters 
and gives rise to non-dimensional parameters of physical significance such as 
the Reynolds number, Grashof number and Prandtl number, etc. which can 
be varied independently to characterise the results. 

The nondimensionalisation of the governing equations and the corre- 
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Fi^re 2.3: An alternate »mpntational domain and boundary conditions 
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spending boundary conditions is achieved by introducing dimensionless in- 
dependent and dependent variables. The definition of these dimensionless 
variables are based on a choice of a characteristic dimension and other char- 
acteristic quantities that depend on the geometry of the physical domain 
and the physical flow and heat transfer properties of the fluid. There is a 
freedom to select these quantities, and the proper nondimensionlisation is 
problem dependent. There may be many possible non-dimensional forms, of 
which some may be more natural and sensible than the others. In this section, 
the dimensionless variables, the dimensionless governing equations and the 
dimensionless boundary conditions, both for natural and mixed convection 
flows, are presented. 

2.5.1 Dimensionless Variables 

The nondimensionlisation of the governing equations and the boundary con- 
ditions is done here by introducing dimensionless variables defined as follows: 


9 = 


X 

X 


^ H 

H 

u 

u 



V= — 

^u'e 

Uc 

T 

-Too 

p_P^-P 

Tyj 

-Too 

pU^c 


where the characteristic velocity, Uc is given for natural convection flows by 


Uc = [9PH{Tw-T^)]K (2.6) 

and for mixed convection flow by 


Uc = Uo, ( 2 . 7 ) 

where Uq is the uniform inflow velocity. 

In the definition of dimensionless variables presented above, the height of 
the channel, H, is employed as the characteristic length dimension, {Tyj—Too) 
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as the characteristic temperature difference and pUl is the characteristic pres- 
sure. The nondimensional velocities are differently interpreted for natural 
and mixed convection according to the different definitions of characteristic 
velocity given in equations (2.6) and (2.7). 


2.5.2 Dimensionless Governing Equations 

With these dimensionless variables, the dimensionless equations of mass, mo- 
mentum and energy transport for natural convection and mixed convection 
can be written separately as follows: 


'Natural Convection: 
mass: 




X-momentum: 


U— + V— = — ^ 

dX dY dX \/GTfj 


d^U &^U 
+ 




Y-momentum: 


dY dP 1 
dx dY dY ^ 


\d^V d'^V 

+ 


dX2 dF2 


energy: 


+ V~ = ^ 

dX dY Pr^/Grff 


dH dH 

+ 


dx2 ar2 


( 2 . 8 ) 


(2.9) 


( 2 . 10 ) 


( 2 . 11 ) 


where Grff and Pr are the Grashof number and Prandtl number, respectively, 
given by 


1/2 


( 2 . 12 ) 
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The corresponding equations for mixed convection are : 
mass; 


dx'^ BY 


(2.14) 


X-momentura: 


du du _dp 1 d'^U' 

^ dY ~ dX'^ Ren dx^^W^ 


Gr 


H , 




Y-momentum: 


,,dV dP 1 

U h V — = 1 

dX dY dY Ren 


'd^v av 


energy; 


de ydd _ 1 r d'^e dH ' 

^dX^^ dY~ PrRen dX^ ^ dY^ 


where Ren is the Reynolds number given by 


Ren = 


UqH 

V 


(2.15) 


(2.16) 


(2.17) 


(2.18) 


The dimensionless equations are seemingly different for natural and mixed 
convection flows. However, although the nondimensionalisation has changed 
the appearance of the equations, the fundamental governing equations are 
the same in both cases. It is merely more appropriate to use different nondi- 
mensionalisation for the natural and mixed convection cases, because of the 
different physics involved. 


2.5.3 Dimensionless Boundary Conditions 

If the aspect ratio of the channel and the dimensionless length of the extended 
domain are, respectively, defined by 



(2.19) 
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= (2-20) 
the dimensionless boundary conditions can be expressed as follows: 


at X = 0, 


.. 1 dU_dV_^ 

dX'^'dX ^^dX 


0,P = 0 


1 cfr 

at ^ = 0)^ 0<X<1: C/ = 0,y = 0,^ = 1,^ = 0 

1 dU ^ ^ 69 ^ dP 


at X — 1 €i, 

M = 0, V = 0, 0 = 0, P = 0 for natural convection 
P = — 1, y = 0, 0 = 0, ^ = 0 for mixed convection 

2.6 Stream Function- Vorticity Formulation 

The dimensionless governing equations in the primitive variables for two- 
dimensional flows have two components of velocity, the pressure and the 
temperature as dependent variables. For incompressible flows, it is a crucial 
fact that there is no direct equation for the pressure — the pressure is merely 
a hnk between the continuity and momentum equations. Primitive variable 
solution procedures usually employ the continuity equation as a constraint 
to derive the correct pressure field through a cumbensome iterative method. 
Keeping this in view, the present problem has been formulated using the 
stream-function vorticity approach, which is presented below. 

The stream function-vorticity formulation is a commonly used approach 
for the numerical simulation of steady two-dimensional or axisymmetric in- 
compressible flows. This formulation has the advantage that it eliminates 
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the pressure from the governing equations, and the stream function automat- 
ically satisfies the continuity equation. Thus, this method has the advantage 
of dealing with only two equations, while the primitive variable formulation 
requires that the continuity equation also be solved along with two momen- 
tum equations. 

In addition, Dirichlet conditions are often given for the stream function, 
instead of the Neumann conditions generally needed for the pressure. On 
the other hand, there is no direct boundary condition for surface vorticity, 
which has to be specified in terms of the stream function field. The accuracy 
of the overall solution critically depends on how this boundary condition is 
handled. In the rest of this section, the stream function - vorticity formula- 
tion is presented along with the boundary conditions. 


2.6.1 Vorticity Transport and Stream Function Equa- 
tions 

The vorticity transport and stream function equations can be obtained by 
introducing the dimensionless stream function, and vorticity, Q, in terms 
of their dimensional forms ip and ui: 


^ = 


jL- 

UcH 


Q = 


uH 

Uc 


These can be related to the non-dimensional velocity field by: 


8Y' 


y = 


ay 


and 


“ dX dY 


( 2 . 21 ) 


( 2 . 22 ) 


We can readily see that Equation (2.21) satisfies the continuity equation 
(2.8) and (2.14). Thus, if we substitute the dependent U,V variables with 
the stream function we need not be concerned with the continuity equation 
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anymore. Invoking equations (2.21) into (2.22) we obtain a Poisson equation 
for the stream function as 


arrays 


(2.23) 


This equation is applicable irrespective of type of flow (natural, mixed or 
forced). 


By differentiating the X and Y momentum equations, (2.9) and (2.10) 
with respect to Y and X, respectively, subtracting the resulting equations 
from each other and rearranging, we obtain the vorticity transport equation 
for natural convective flows as 


dO ydO _ 1 fd^Q d‘^n\ 89 

^ax ''' ''ay ~ \sx‘> ay^)'*' ar 


(2.24) 


Similarly, the vorticity transport equation for mixed convective flows can 
be obtained as 


80 80 ^ 1 (cf^O 8‘^0\ Gm 89 

8X'^ 8Y Ren [dX^ 8Yy Re^ 8Y 


(2.25) 


2.6.2 The Pseudo-Trsinsient Forms 

Instead of directly solving the stream function, vorticity and energy equations 
(2.23, 2.24 or 2.25, and 2.11 or 2.17), we solve the pseudo-transient forms 
of these equations as a solution strategy. The stream-function equation is 
written in the Cauchy-Kowaleska form : 


8t 


92^ 32^ 


(2.26) 


which is solved along with the following equations: 
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for natural convective flow; 

dn dUQ dvo _ 1 r , ae 

dr'^ dx ^ dY ~ [dx'^ ^ ^ 


de due dve _ i ( an \ 
dr^ ax ^ aY ~ Pry/GFF ’^dY^J 

or, for mixed convective flows: 


ao auQ avo _ i a^f^l Orn de 

dr ^ ax av ~ Res [dX^ ^ aF^J + Rejj dY 


(2.28) 


(2.29) 


ae due dve _ i r a‘^e , d'^e 1 

ar ax ^ ay ~ RenPr [ax2 ^ ay^J ^ ^ 

Note that the governing equations have now been put in the conservative 
form. In these equations r is a pseudo-time-like variable. The equations 
are solved to obtain ‘unsteady’ solutions which converge to the final steady- 
state. The transient terms involving r reduce to zero at the steady state, 
leaving solutions of the original steady-state equations. The use of these 
pseudo-transient forms of the equations allows us to obtain the steady state 
solutions easily and naturally by time-stepping from arbitrary initial condi- 
tions. Other iterative methods may be more efficient in obtaining the steady 
state solution, but would probably be more uncertain, given that three cou- 
pled partial differential equations are involved. Furthermore, tailoring these 
other iterative schemes to the special, full-matrix, spectral element formula- 
tion is likely to prove cumbersome. 


2.6.3 Boundary Conditions 

The vorticity transport and the stream function equations require bound- 
ary conditions to be specified on all the boundaries of the solution domain. 
The boundary conditions for the primitive variables formulation have already 
been discussed in Section 2.4.2. In this section, the boundary conditions for 



46 


CHAPTER 2 . MATHEMATICAL FORMULATION 





2.6. STREAM FUNCTION-VORTICITY FORMULATION 


47 


the stream function and vorticity are presented. Figure 2.4 depicts the com- 
putational domain, with the boundary conditions used. 

The impermeability condition at the wall AB implies that there is no fluid 
flow across this boundary and that the solid boundary itself can be considered 
a stream line. Obviously, the stream function is constant along this bound- 
ary, and as one value of the stream-function can be assigned arbitrarily, ^ = 
0 is conveniently specified along this boundary. It also follows that on this 

boundary = 0. This can be used to obtain the boundary condition of 

vorticity (see, e.g.. Comini et al., 1997) by substituting into Equation (2.23) 
^ 2 ^ 

to get D, = — at the solid boundary AB. A similar boundary condition 
for vorticity is used on the other wall EF. 

The boundary EF corresponds to a solid boundary and, therefore, the 
stream function is constant along this boundary due to the reason already 
mentioned above. Therefore, a general condition ^ is specified along 
this boundary. The value of is intimately related to the upward volume 
flow rate per unit depth, v, in the channel which is equal to the difference 
(= 0 — ■^ft) between the stream function values at the two plates. For 
mixed convection flows, the volume flow rate is imposed and known, and the 
value of 'F(i can be specified in advance. However, for natural convection 
flows, the volume flow rate is induced and not known a priori, and has to 
be determined iteratively as part of the numerical solution procedure. 

The symmetry across BC allows no flow to cross this boundary which, 
therefore, is a stream line. It is clear from Figure 2.4 that the stream lines 
corresponding to boundaries AB and BC are the same. Thus, the value of 
the stream function along the boundary BC is also specified as zero. The 
value of the vorticity Q. is equal to zero on this boundary because of symme- 
try. The symmetry boundary DE is similar to BC except that the stream 
function value is tfie same as on EF. 

At the inlet boundary CD, the flow is assumed to enter at the ambient 
temperature (^=0) and with zero vorticity, as the flow is as yet away from 
the walls which create vorticity. The stream-function boundary condition 
is obtained by assuming that the inflow is vertical (i.e., 1^=0), yielding the 
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Neumann condition || =0 along this boundary. 

At the outflow boundary AF it is again assumed that the flow is vertical, 
go H =0. It is also assumed that the flow is fully-developed and the nor- 
mal gradient of the temperature and vorticity are thus taken as zero. These 
boundary conditions frequently appear in the literature. Even though fully 
developed flow may not be achieved at the exit of the channel, the zero gradi- 
ent boundary conditions generally allow the flow distribution to be obtained 
quite accurately (Huang and Vafai 1994). 

To implement the vorticity boundary condition at the solid walls we do 
the following. Using a Taylor series expansion of the stream function on 
the solid boundary and incorporating the additional condition = 0, the 
vorticity at the solid boundary is evaluated by 


— 


Ay2 


at y = 0, 1 


where Ay is the grid spacing at the solid boundary. The stream function 
value at the wall is 0 or on the left and right walls, respectively, and 
is the value at the corresponding neighbouring interior point. 


Mass Flow for Natural Convection 

As mentioned above, the volume flow rate, hence ^6, in natural convective 
flows is not known in advance, and needs to emerge as part of the solution. 
This can be done by ensuring that the momentum change along the vertical 
channel is exactly balanced by the combined viscous shear forces at the walls 
and the buoyancy forces on the fluid in the channel. A straight-forward 
application of the integral momentum equation at the inlet and outlet of the 
channel jdelds the following equation for the pressure drop; 



(2.31) 
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where is the difference of the average pressure at the inflow and outflow 
boundaries, which should be zero for natural convection. 

The flow rate, and hence clearly depends on the Grashof number, 
Grfr, and the aspect ratio, Ar, through an unknown relation. For the cor- 
rect value of Wb, buoyancy and shear forces will be exactly balanced by the 
momentum change in the channel, and AP^t, will be zero. In the iterative 
procedure for natural convection, the value of is guessed and the steady 
state solution is obtained. By using the steady state solution in Equation 
(2.31) it is determined from if AP^^ = 0. If not, is changed accordingly 
and a new steady state solution is obtained, and so on till APo„ = 0 (to six 
decimal places). 

It is clear that the natural convection case is very computer-intensive, as 
it involves obtaining many steady state solutions for each case, while mixed 
convection requires only one solution for each case. 


2.7 Criteria for Optimum Spacing 

In this section, the use of the numerical solutions in obtaining the optimal 
spacing (for maximum heat transfer) for the multi-channel problem discussed 
in Section 2.2.1 is elaborated upon. 

For multiple channels formed by symmetrically heated equal-spaced par- 
allel plates maintained at the same constant temperature, the flow and heat 
transfer characteristics of each channel would be the same. Therefore, it is 
sufficient to analyse only one channel to obtain results pertinent to a series of 
identical vertical parallel-plate channels. Prom the numerical solution of the 
problem depicted in Figure 2.4, whether for natural or for mixed convection, 
we obtain the steady state velocity and temperature fields in the solution 
domain. Prom the temperature field we can obtain the heat transfer from 
the plates in the following way. 

We define the heat transfer coefficient at the walls with reference to the 
temperature difference (T^-Too). The local heat transfer coefficient is defined 
as: 
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(2.33) 


where the local heat flux at the wall, Qv,, can be obtained from the tempera- 
ture field by 

<lw = -k 

where k is the thermal conductivity of the fluid. 

By non-dimensionalising this equation, we obtain 

(T.-Too) fde\ 




Qw 


= -k- 


H 


from which the local wall Nusselt number, 
Nu = — = 


.dY, 


'89' 

,ay 


(2.34) 


(2.35) 


can be obtained, using the numerical solution for 9. 

The average wall Nusselt number can then be obtained as 

1 U 

I^Uav = T? / Nu dx— Nu dX, 

H Jo Jo 


(2.36) 


which can again be evaluated from the numerical solution. The total heat 
transfer (per unit depth) can then be obtained by 


Qw — havH (Tyj — Too) — kNuav {Ty, — Tqo). 

where hav is the average heat transfer coefficient. 


(2.37) 


When there are multiple channels, the number of such channels that can 

be fitted into a span W is We thus see that the total heat transfer in 
that span will be 

^ W Wk 

Q,.m = 2 - = 2 - T „) 


(2.38) 
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where Ar {= H/L) is the aspect ratio of each channel, and the factor of 2. 
arises because there are two heat transfer surfaces per channel. 

We now see that, for fixed W, H, (T^ — Too), and k, 


Q total OC Nu av^r (2.39) 

which provides us the objective function (=A^Ua„Ar) for our optimisation 
problem of finding the maximum heat transfer for the case of multiple chan- 
nels. 

Thus, for fixed Grashof number Grn (for natural convection) or fixed 
Gth and Ren (for mixed convection), the numerical solutions are computed 
for a variety of aspect ratios to determine the value of Ar that yields the 
maximum value of NuavAr- This aspect ratio is then taken to correspond to 
the optimum spacing between the plates, for the multiple channel case. In 
contrast, for the single channel case the optimum Ar is simply the one that 
gives the highest Nuav, which also maximises the heat-transfer for that single 
channel. The interpretation of the numerical results is thus quite simple. 


^ ''Tit 


1 . 450 . 42 . 
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Chapter 3 

Numerical Formulation 


3.1 Introduction 

The fundamental laws of fluid flow and heat transfer are expressed as a set of 
coupled, non-linear partial differential equations, often with complex initial 
and boundary conditions. In most situations, analytical solutions of these 
equations either do not exist or are mathematically complex in form. As 
a consequence, one generally resorts to numerical solutions using a suitable 
discretization method. The discretization methods that are most commonly 
used are the flnite difference method, the flnite element method, and, more 
recently, the spectral method. 

In the present thesis, the numerical solutions of the governing equations 
are obtained by the spectral element method of Patera (1984) implemented 
on the stream-function vorticity formulation by Prakash (1995). The spec- 
tral element method is basically a synthesis of the spectral and the flnite 
element methods (Korczak and Patera, 1986). In this method of discretiza- 
tion, the computational domain is decomposed into a number of elements, 
and within each element the dependent variable is represented as a high-order 
Lagrangian interpolant in terms of Chebyshev polynomials, the coefficients 
of which are related to the function values at the Gauss-Lobatto Chebyshev 
collocation points. The governing equations are partially discretized using a 
flnite difference scheme for the time derivative term leading to a second-order 
elliptic Helmholtz equation. While the convection and the source terms are 
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treated explicitly, the diffusion terms are treated implicitly. Xhe Helmholtz 
equation is further discretized by spectral element method, using the vari- 
ational principle, resulting in a system of linear algebraic equations. These 
systems of equations is efficiently solved by the LU decomposition method. 

In this chapter, the discretization of a typical computational domain, 
and the temporal and spatial discretization of the conservative form of the 
governing equations are presented. The formation of elemental and global 
equations and the implementation of boundary condition are described in 
detail. The solution method of the system of algebraic equations and other 
relevant computational details are provided. 


3.2 Discretization of the Domain 

The first step towards the numerical solution of the mathematical model is 
to discretize the computational domain. In the spectral element method, 
the computational domain is divided into a number of rectangular elements 
spanning the X and Y directions. Each element of the computational do- 
main is further subdivided into a set of grid points. Figure 3.1 illustrates a 
nominal 4-element discretized computational domain. Figure 3.2(a) depicts 
the t 5 rpical element of the computational domain showing {X, Y) coordi- 
nates of its comer points. 

For the sake of convenience in using Chebyshev interpolating functions, a 
local coordinate system applicable to each individual element is introduced. 
Consequently, each rectangular element in the global coordinates (X, T) is 
mapped onto a square element in the local coordinates {x, y) spanning -1 
to +1 in the x and y directions, as shown in Figure 3.2b. Thus, the local 
i element coordinate system in terms of global coordinates system can be 
defined as 


= (3.1) 

t = -^(Y - ai) - 1 

y 


(3.2) 
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Figure 3.1; A nominal discretized computational domain 
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Figure 3 . 2 : A typical element in (a) the original coordinates, and (b) trans- 
formed local coordinates. 

where (aj., a^) and 6 * ) are the coordinates of the top left and the bottom 
right comers of the element, and Lj.(= 6 * 3 . - o^) and Ly{= are the 

lengths of the element in the X and Y directions respectively. It is to be 
noted that the global coordinates {X, Y) used in the mathematical formula- 
tion of the problem can be easily retrieved from the local elemental coordinate 
just by inverting the linear transformation given by equations ( 3 . 1 ) and ( 3 . 2 ). 

Fi^re 3.3 shows the element in the local coordinate system (x, y) hav- 
ing Nl+1 and N^+l number of grid points in x and y directions, respectively. 
The position of any grid point within the element is identified by a set of two 
indices, say (j, k), the top left comer of the element being ( 0 , 0 ), and the 
bottom right comer (iV*, Aj). The local coordinates {x), f^) of the (j, k) grid 
point are generated by the Gauss-Lobatto Chebyshev formula as follows: 


x) = cos(^) 

) — 0 , 1 , 2 , .. 

. 5 i 

(3.3) 

-i 

Vk = cos{—) 

fc = 0 , 1 , 2 ,.. 


(3.4) 
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( 0 . 0 ) 



Figure 3.3; The grid-points and local two-index numbering in a typical ele- 
ment. 

3.2.1 Grid Point Numbering 

Local Numbering 

For the sake of convenience in writing the discretized elemental equations in 
the conventional matrix form, the grid points identified by two indices, say 
(j,k), are also specified by a single index, say p, using the transformation 
given by 


p = (iVy -f l).j + k i = 0, 1, ..., iV^, k = 0,l, ■■■, Ny (3.5) 

where, and Ny+1 are the number of grid points in X and Y directions 

respectively. 

Figure 3.4a illustrates a typical element of the computational domain hav- 
ing 4 grid points in the X direction and 5 grid points in the Y direction: in 
this figure, every grid point is identified by two indices, i.e., (y, k). Beginning 
from (0, 0) at the top left corner, the successive neighbouring grid point is 
marked by incrementing either j or by 1 in accordance with its location 
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Figure 3.4: The (a) two-index and (b) one-index local numbering of grid- 
points in a t 3 T)ical element. 


with respect to its neighbour. The corresponding single index numbering of 
the grid points at element level is shown in Figure 3.4b. Starting from 0 at 
the top left corner of the element, the grid points at the element level are 
numbered in serial order along the Y direction by advancing line by line in 
the X direction. 


Global Numbering 

The numbering systems above are local to each element, as every element will 
have the same set of numbers [(0, 0),..{N^, Ny) or 0, 1, ..., {Nj;+l){Ny 4- 1) - 1] 
describing its grid-points. Apart from the local two- and one-index number- 
ing, the domain also has a global numbering system in which each grid point 
has a unique number across all elements. This facilitates the formation of a 
global equation set by combining the elemental equations. Thus every grid 
point has two numbering systems associated with it. The local numbering 
corresponds to its position within an element, and the global numbering to 
its position in the computational domain as a whole. 
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Figure 3.5: The global numbering of grid points in the computational domain. 


Figure 3.5 depicts the global grid point numbering system of the compu- 
tational domain consisting of 2 elements in the X direction and 2 elements 
in the Y direction. Each element of the computational domain is assumed 
to have 4 grid points in the X direction and 5 grid points in the Y direction. 
The elements are identified by two indices, say {ie,je), marked at the centre 
of each element and the grid points by single global number. Starting from 
1 at the top left corner of the computational domain, the grid points at the 
global level are numbered serially along the Y direction by advancing line by 
line along the X direction. 
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3.3 Discretization of Governing Equations 

Now that we have discretized the computational domain, the next step in- 
volves the equation discretization which approximates the governing differen- 
tial equation by a system of linear/nonlinear algebraic equations in variables 
defined at the grid points of the computational domain. For transient/ pseudo- 
transient problems, the equation discretization involves two components, the 
temporal discretization and the spatial discretization. In the succeeding two 
subsections the temporal and spatial discretization of the conservative form 
of the governing equations are presented. 


3.3.1 Temporal Discretization 

In the spectral element method, the temporal derivative term is generally 
treated by a finite difference scheme. Both explicit as well as implicit schemes 
have been exploited in a variety of heat and fluid flow situations. The ma- 
jor disadvantage of the explicit schemes is the restriction imposed on the 
maximum time step by the numerical stability criterion. In contrast, the 
implicit schemes are usually free from these restrictions and, therefore, allow 
relatively larger time steps as compared to explicit schemes. However, with 
regards to the number of arithmetic operations required per time step, the 
implicit schemes are typically more expensive. Implicit formulations also of- 
ten have greater algorithmic complexity. 

Although, explicit schemes are only conditionally stable, this is not al- 
ways bad. In fact, the standard stability criterion for explicit schemes applied 
to convective equations, the CFL (Courant-Priedrichs-Lewy) condition, con- 
veys a clear physical meaning. It states that the time step should not be 
greater than the time taken for any fluid particle to cover a distance equal 
to one grid interval. It is observed that even for implicit schemes that are 
unconditionally stable the accuracy of numerical solution deteriorates if the 
CFL condition is not satisfied. Obviously, it is always safer to satisfy the 
CFL condition for strongly convective flows and, therefore, it is unnecessary 
to implement complex implicit schemes for the convective terms. On the 
contrary, for the diffusion terms it is usually more efficient to implement im- 
plicit schemes because of their better stability properties. 
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In situations when both convection and diffusion terms are present, it is 
often the practice to use an explicit representation for the convective terms 
and an implicit representation for the diffusion terms. For example, in the 
Adams-Bashforth-Crank-Nicolson scheme, the convective terms are treated 
by the explicit Adams-Bashforth scheme and the diffusion terms by the im- 
plicit Crank-Nicolson scheme. Both of these schemes are second order accu- 
rate, and the overall scheme also possesses second order accuracy in time. 
However, it is only conditionally stable as well as algorithmically involved. 

Patera (1984) and Korczak and Patera (1986) adapted the Adams-Bashforth- 
Crank-Nicolson scheme for a spectral element temporal discretization. In 
their formulation, the implicit representation for the diffusion term is the 
keystone of the method — as it allows the basic governing equations to be 
expressed in terms of a Helmholtz equation which can be solved using the 
variational principle. Their scheme is second order accurate in time but is 
exponentially accurate in space because of its spectral spatial discretization. 


In the present thesis, the pseudo-transient forms of the governing equa- 
tions are discretized in time by treating the diffusion term implicitly and 
remaining terms explicitly. Although, the temporal discretization used is 
only first order accurate, this does not affect the steady state solution of the 
problems in which we are actually inter^ted. The time discretized form of 
the general transport equation 


84) , (d4>U d4>V\ 


= r 


d'^4> 


+ S,p 


is given by; 


^n+l _ 

At 


+ A 


d4)U d4)V 


dx dY 


= r 


8^ ^ 8‘^(l> 


8X^ 8Y'^ 


n+l 


+ 5? 


(3.6) 


Equation (3.6) can be rearranged to obtain the following form: 


8X'^ j ^ PAr 


4)^ A \84)U 84)V 

PAr ^ r 8X 8Y 


5 'n 
4> 

T 


(3.7) 
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or 


(l>xx + <^VT “ — / 


(3.8) 


Equation (3.8) is the well-known Helmholtz equation. The function (j) is 
the unknown to be computed and / is the known forcing function. In equa- 
tion (3.6) and (3.7), the superscript n refers to the known values at current 
time r, and n+1 to the unknown values at time r + At. The parameter 


A^(= -4-) can take different positive values. In the limiting case when 
FAr 

0, equation (3.8) reduces to the Poisson equation which too can be 
solved using the Helmholtz solver. 


3.3.2 Spatial Discretization 

In the spectral element method, the solution of the Helmholtz equation is 
based on the variational formulation requiring the minimization of a corre- 
sponding functional. 


3.3.3 Variational Formulation 

The functional corresponding to the Helmholtz equation [Equation 3.8] sat- 
isfying homogeneous Neumann boundary conditions can be obtained as 




dXdY 


(3.9) 


where the integration spans over the entire computational domain. 


When the computational domain is decomposed into M rectangular ele- 
ments, the functional I can be simply expressed as the sum of the functional 
of the individual elements, ie., 


M 


i=l 


(3.10) 


where I is the functional corresponding to element and can be expressed 
as 
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// 


+ (4-r) - 


dXdV 


(3.11) 


the integration spanning over the area of the i*'' element, defined as the rect- 
angle between [ 0 ^, 0 ^] and [bl.,by]. 


If there are {N^+l) (Ny+1) Gauss-Lobatto Chebyshev grid points within 
the element, the dependent variable <j>{X, Y) and the forcing function, 
f{X,Y), can be expressed as Lagrangian interpolants in Chebyshev polyno- 
mials of the local coordinates,—! < 5*,^ < H-l, as 




j=0 k=0 


(3.12) 


and 


Ni Ni 

5‘) = E E (3-13) 

j=0 fe=0 

where and /j^. are the values of <() and / at the grid point (x}, yl) respec- 
tively. The interpolation functions hlj{x^) and h\{y'') are Nl—th and N^ — th 
order local Lagrangian interpolants satisfying 


hj{xl) = 5ji 

and 


Kiyln) = ^km 


(3.14) 


(3.15) 


at the local grid points within the element, where Sji and Skm are the 
Kronecker delta functions, xj is the coordinate of Ith. grid point in the x*- 
direction and ^ is the coordinate of the mth grid point in the y^-direction. 


The interpolation functions /i}(x*) and hUf) can be expressed in terms 
of Chebyshev polynomials as 
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or 


4>xx + <^yv ~ — f (3.8) 

Equation (3.8) is the well-known Helmholtz equation. The function (j) is 
the unknown to be computed and / is the known forcing function. In equa- 
tion (3.6) and (3.7), the superscript n refers to the known values at current 
time r, and n+1 to the unknown values at time r + At. The parameter 

A^(= p^) can take different positive values. In the limiting case when 

0, equation (3.8) reduces to the Poisson equation which too can be 
solved using the Helmholtz solver. 

3.3.2 Spatial Discretization 

In the spectral element method, the solution of the Helmholtz equation is 
based on the variational formulation requiring the minimization of a corre- 
sponding functional. 

3.3.3 Variational Formulation 

The functional corresponding to the Helmholtz equation [Equation 3.8] sat- 
isfying homogeneous Neumann boundary conditions can be obtained as 


r = 1 1 -jCWx)" + (•hf} - H 


dXdY 


(3.9) 


where the integration spans over the entire computational domain. 


When the computational domain is decomposed into M rectangular ele- 
ments, the functional I can be simply expressed as the sum of the functional 
of the individual elements, ie.. 


M 

i = Y.r (3.10) 

i=l 

where P is the functional corresponding to element and can be expressed 
as 
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/' = // [-|{(4)" + (4)^ - )“ - 4‘‘A dXdY (3.11) 

the integration spanning over the area of the element, defined as the rect- 
angle between [a^ja^] and [bl,bl]. 

If there are (A^*-f-l) (Ar*-|-1) Gauss-Lobatto Chebyshev grid points within 
the element, the dependent variable <j){X,Y) and the forcing function, 
f{X,Y), can be expressed as Lagrangian interpolants in Chebyshev polyno- 
mials of the local coordinates,—! < < -t-1, as 

Ni Ni 

^ (3.12) 

j=0 fc=0 

and 

(3.13) 

j-Q k=0 

where and /jj^. are the values of (f> and / at the grid point (x), y^) respec- 
tively. The interpolation functions /i}(x‘) and hKy^) are N].—th and Ny — th 
order local Lagrangian interpolants satisfying 

hij{xl)=5ji (3.14) 

and 

(3.15) 

at the local grid points within the element, where dji and 5km are the 
Kronecker delta functions, x\ is the coordinate of Ith. grid point in the x'- 
direction and is the coordinate of the mth grid point in the y*-direction. 

The interpolation functions h)(x‘) and can be expressed in terms 

of Chebyshev polynomials as 
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2 1 


and 


71=0 


N* 

>‘i(s‘) = 4 E ;^r„(i/i)r„(r) 


m=0 ^kCm 

where r„ are the order Chebyshev polynomials given by 


(3.16) 


(3.17) 


Using Equation (3.3), 


and 


Tn(x^) = cos(ncos 

(3.18) 

T^(^})=cos(^) 

(3.19) 

Tmivl) = COS{ ) 

(3.20) 


The C's in equations (3.16) and (3.17) are the integer functions given by 


and 


Cj = 1, 

for j ^0 or Nl 

(3.21) 

= 2, 

] 

0 

11 

0 

(3.22) 

1 

Ck = l, 

for k^O or Ni 

(3.23) 

= 2, 

II 

0 

0 

<cj 

(3.24) 

Now, the equations (3. 

12) and (3.13) can be rewritten as 



^ Ni Ni AT* 


S‘) = E E E E (3.25) 

^ '‘x-^ '‘y jzzO k=<} n=0 m=0 

TUs'Mf) 
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and 


Ni 


^ x: £ E E 4^ 7^ Ux^)Tn{x^) (3.26) 
J V / j^Qi.=0n=0m=0 ^ 

Tmiyl)Tm{.f) 


As the physical rectangular elements are mapped into square elements 
using the linear transformation given by equations (3.1) and (3.2), the dit 
ferentiation or integration performed in the transformed coordinates (5%^^) 
must be properly scaled. As a result, the corresponding derivatives in ttr^ 
two coordinate systems are related by 


and 


dp _ 2 dp 

Jx ~ Li dx^ 

(3.27) 

dp 2 dp 

W^^UyW 

(3-28) 


where, ^ and ^ are the scale factors. Therefore, the derivative terp^g 


Li 


consisting of (j)^x and (py iu equation (3.11) can be easily obtained as 
. , Ni^i 1 

* = 4 E r S S 'l>)krrrr^M^i}TM) 

^ LlN^N^ ^Q^Qn=Om=0 CjCnCkCm 

dxi 

and 


(3.29) 


Ni Ni 


4 = 


LlNiNl^,^o 


r, T^n Jfc=0 n=0 m=0 OjUnOkUm 

^ ,_i.dTrr,{f) 


(3-30) 
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3.3.4 Formation of Elemental Equations 

The elemental algebraic system of equations are obtained by minimizing the 
functional given by equation (3.11) which requires that the variation of P 
with respect to the nodal values (pjj. vanish, i.e., 


5P 




= 0 


A 


{3.31) 


At this point, it is important to note that it follows from Equation (3.25) and 
the following two discrete orthogonality properties of Chebyshev polynomials 
on a Gauss-Labatto grid: 


Ni 


- y r„(x')r„(x‘) = c^CA 


n=0 


PJ 


and 

TV* 

m=0 

that 


d 

~ (3.32) 

Inserting the 4>^ and p jfrom Equations (3.25) and (3.26) and the deriva- 
tives (j)\ and (py from Equations (3.29) and (3.30) into the functional given 
by Equation (3.11), performing the resulting integrals and finally applying 
the Equations (3.31) and (3.32), the elemental algebraic system of equations 
will simply be (henceforth the Einsteinian repeated index convection will be 
used): 
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Here. 


4 ' 


Li 


Li 


1 7* 

t 

i\2 A km JJ 


:Al^m 


fi fi 1 

R* = Z1—JL———RX UV 
^jklm ^ A-^km. 


Ajl = P P ^ ^ ^ p Tn{x j)Tm.{Xi)o.r, 

n=0m=0 k^rik.'m 


B 


A ^ ^ 1 


EE 




-Tn{x^^TUxi)br, 


and 


- L 


1 dTndTm^-i 


-1 dx} dx^ 

0, n + m. odd 

~~^[J\{n—m.)/2\ ‘^|(n4-m)/2|]i d” GVBn 


where. 


and 


= 0 , k = Q 


bnm = j ^ TnTmdx' 


= 0 . 


n + m odd 


1 1 
+ 


1 — (n + m)'^ 1 — (n — m,)^ ’ 


( 3 . 35 ) 

( 3 . 36 ) 

( 3 . 37 ) 

( 3 . 38 ) 

( 3 . 39 ) 

( 3 . 40 ) 

( 3 . 41 ) 

( 3 . 42 ) 


n + m even 
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Note that the matrices and B]i are independent of the elemental 
length Li, L\ and are the same for all the elements, which have the same N* 
and Ni- Thus, these universal matrices and can be computed first 
and then the individual elemental matrices and can be computed 
for every individual element. 

Equation (3.33), with four-dimensional arrays B and C, is not convenient 
for direct numerical solution. Consequently, this elemental equation is trans- 
formed into more conventional form by combining the two indices j and k to 
give a single index p, as in Equation (3.5). Similarly, the other two indices I 
and m are combined to give another single index q. As & result, the equation 
(3.33) is transformed to conventional matrix form given by 

= Byi (= ej,) (3.43) 

For the element having {Nx+l){Ny + l) grid points, the coefficient matrix 
[C] will be of the order of {Nx + l){Ny 4- 1) x {Nx + l){Ny -f- 1). 

3.3.5 Formation of Global Equations 

In the preceding section, we obtained the system of discrete equations for each 
element. Once the elemental equations for all the elements of the computa- 
tional domain are obtained, these are combined to form the global system 
of equations by adding the stiffness contributions corresponding to common 
boundary nodes of the adjoining elements. The formation of global matrix, 
[Cgi°^], from elemental matrices, [C], and right side global vector, 
from elemental vectors, {e}, are accomphshed as follows: 

For every grid point that lies exclusively within an element, the coefficient 
of the elemental matrix [C] in a locally numbered row and column are trans- 
ferred to the corresponding globally numbered row and column in the global 
matrix it ig the other hand, the coefficient of the global ma- 

trix corresponding to a grid point at an elemental interface is obtained by 
adding the elemental coefficients corresponding to that grid point from all the 
elemental sets in which it appears. However, in doing so, the correspondence 
between the local and global number of the grid point in the computational 
domain is maintained. The same procedure is adopted for the formation of 
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The total number of grid points in the computational domain consisting 
of r.T: elements in the x direction and Vy elements in the y direction, each 
element having N^ and Ny grid in x and y directions respectively, will be 
i Nxfx + + 1) corresponding global coefficient matrix 

will be of the order (N^Vx + l){Nyry + 1) x {NxTx + l){Nyry + 1). The global 
equations then are 


Qglobal global _ ^global (3 44 ') 

where, the vector incorporates all the unknown (pim of all the elements 

and the indices pg and Qg represent the global indices given by 

Pg = l,...,{Nx.rx + l).iNy.ry + l) (3.45) 

Qg — 1 , •••, {Nx-Tx + 1) .(^Ny .Ty + 1 ) 

3.4 Imposition of Boundary conditions 

The boundary conditions corresponding to the stream function- vorticity for- 
mulation of the problem presented in section 2.7.2 are basically of the two 
t\-pes; homogeneous Neumann boundary conditions and Dirichlet boundary 
conditions. The variational formulation of the Helmholtz equation (3.8) using 
the functional (3.9) automatically satisfies homogeneous Neumann boundary 
conditions and, therefore, no modification of the elemental coefficient matri- 
ces in Equation (3.3.3) is needed to impose these boundary conditions. 

Dirichlet boundary conditions may be imposed by matrix condensation. 
Accordingly, the rows and columns corresponding to grid points ha^dng Dirich- 
let boundary conditions axe eliminated from the system matrix. Although, 
this method leads to computational economy by reducing the order of the 
original system matrix, it is somewhat complicated to implement. 

To impose Dirichlet boundary conditions we apply a simple stratagem 
commonly used in the finite element method. The diagonal element of global 
matrix in Equation (3.44) corresponding to the grid points having 

Dirichlet boundary conditions are equated to a very large number ( 10^^°, say) 
nid the corresponding element in the right side vector, are set equal 
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to the specified boundary value multiplied by the same very large number. 
Then the resulting global equations are solved to obtain the solution. The 
advantage of this method lies in its simplicity of implementation, at the cost 
of only slight inefficiency of computation. 


3.5 Solution of Global Equations 

Once the global system of algebraic equations (3.44) is suitably modified to 
impose Dirichlet boundary conditions, an efficient method for solving linear 
system of algebraic equations is needed. Generally, the selection of an effi- 
cient matrix solver depends on the structure of matrix equation. 

In the present case, the left hand side coefficient matrix is typically sparse. 
In addition, only the right hand side vector in the global matrix equation 
takes on different values in the successive time levels. For this situation, 
where the left hand side coefficient matrix remains the same in each succes- 
sive time level, LU decomposition has a distinct advantage over other matrix 
solution methods. For, once we decompose the coefficient matrix into lower 
(L) and upper (U) triangular matrices, we can solve the matrix equation 
(3.44) for different right hand side vectors as many times as required to ob- 
tain the state solution merely at the cost of one forward and one backward 
substitution for each time step. Thus, the LU decomposition, being compu- 
tationally least expensive, is employed to solve the system of linear algebraic 
equations (3.44). 


3.6 Numerical Solution Procedure 

The temporal discretization of the general transport equation for the flow 
property (j) has been presented in section 3.3.1. Adopting exactly the same 
discretization scheme, the time discretized form of each of the governing 
equations for temperature, vorticity and stream function are represented by 
a Helmholtz equation as given in Equation (3.8). The corresponding param- 
eter and forcing function / are defined in Tables 3.1 and 3.2 for natural 
and mixed convective flows respectively. 
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The spatial discretization of the Helmholtz equation (3.8) by spectral el- 
ement method leading to a system of linear algebraic equations (3.44) has 
been presented in section 3.2.2. Each of the Helmholtz equations for tempera- 
ture, vorticity and stream function is discretized to obtain the corresponding 
system of linear algebraic equations. As the governing equations for tem- 
perature, vorticity and stream function are coupled, these systems of linear 
algebraic equations are solved in tandem at each time step. The asymp- 
totic steady state solutions of the governing equations are obtained by time 
marching. The procedure is as given below: 

1. Specify the initial guess (n=0) for the temperature, vorticity, stream 
function and velocity components. 

2. Solve the energy equation (2.28 or 2.30) for the (n-f l)th time step. 

3. Solve the vorticity transport equation (2.27 or 2.29) for the (n-l-l)th 
time step. 

4. Solve the Cauchy-Kowaleska equation (2.26) for the stream function 
for the (n-l-l)th time step. 


5. 


Compute the velocity field from the stream function using Equation 

( 2 . 21 ). 


6. Check the convergence by the criterion 


rms = 




^total 


< e 


for all fields. 


(3.46) 


7. Update all dependent variables, i.e., *=1: ••• Ntotai 

Step 2-7 are repeated until the condition (3.46) is satisfied by all dependent 
variables (f>. 


3.7 Spatial Derivatives at Collocation Points 

Like the spectral method, spectral element method also allows an accurate 
evaluation of derivatives of dependent variables. The spatial derivatives of 
the dependent variable (j) with respect to X and Y are given by equation 
(3.29) and (3.30). For any point {X = Xr, Y — in the element, 
equation (3.29) can be rewritten as 



O O ^ X 

rw^T.'i'ipFniDP, 

j=0n=0 


(3.47) 
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Table 3.1 Expression for Parameter and forcing function / 
(natural convection flow) 


Equation 

4> 


1 1 

Energy 

9 




Vorticity transport 

Q 



W- 

Stream function 


1 

At 



Table 3.2 Expression for Parameter and forcing function / 
(mixed convection flow) 


Equation 

4> 

A' 

7 ^ - 

Energy 

9 

PtRen 

At 

~PrRe„ £-(g + ^ 

Vorticity transport 

Q 

fieg 

At 

_ Rf>„ - I ^ 4 . svn Gtu d 9 

At - Rif-gx j 

Stream function 





where 


and 


Fnj 


1 


C^Cn 


■^cos 



(3.48) 


= 


(^) 

sin 


te) 


r ^ 0 or Nx 


(3.49) 


Similar treatment works to evaluate the derivative with respect to Y. 


3.8 Accuracy of the Simulations 

In the simulations that follow in Chapter 4, the computational domain is 
t)T)ically discretized in 5 x 3 or 6 x 2 (and for hig h aspect ratios, 6 x 1) 
elements. The number of grid points in each element was always iVj;=7 and 
the two directions. While the number of elements used may not 
seem to be very many, it should be recalled that the spectral element method 
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Table 3.3 Grid Convergence test based on average Nusselt numbers for three 
different aspect ratio of the channel in the natural convection case at Gth = 
10 ® 


A. 

NUav 


5x3 

elements 

6x2 

elements 

5x3 

elements 

6x2 

elements 

1.43 

11.343 

11.395 

11.215 

11.278 

1.52 

11.337 

11.366 

11.208 

11.257 

1.67 

11.311 

11.323 

11.189 

11.240 


has exponential order accuracy, which means that the truncation error de- 
creases faster than any finite power of the number of grid-points. 

To illustrate the accuracy of the simulations, the computed average Nus- 
selt numbers for three different aspect ratios in the natural convection prob- 
lem, at Gth = 10®, are shown in Table 3.3. The Nusselt numbers are com- 
puted in two different ways: Nuav is obtained from the temperature deriva- 
tive at the walls, and from an integral energy balance. The Table 3.3 
shows that the computed values of the two Nusselt numbers are very close, as 
are the results from the 5x3 and 6x2 element discretizations. Throughout 
this study, the difference in the values of Nuav and has been used as 
diagnostic of the numerical error. This difference is typically within 1.5%. 
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NUMERICAL FORMULATION 



Chapter 4 

Results and Discussion 


In this chapter we present the results for the optimum spacing problem 
for the cases of (a) natural convection, and (b) mixed convection, between 
a series of heated vertical parallel plates. The formulation of these prob- 
lems has been already presented in Chapter 2 of this thesis and will be only 
briefly summarized here. The numerical technique which applies the spec- 
tral element method to stream function, vorticity and energy equations had 
already been described in Chapter 3. Below we first present the results for 
the natural convection case, and then consider the mixed convection problem. 


4.1 Natural Convection Case 

The physical problem has been depicted in Figure 2.1 which shows a verti- 
cal parallel plate channel with symmetrically heated walls simulating either 
one element of a parallel array of such channels, or else as a single channel 
standing alone. 

The equations solved for this case are (2.26), (2.27), and (2.28) presented 
in Chapter 2. The false transient technique is used to arrive at the steady 
state solutions of the stream function, '5^, the vorticity, and the temper- 
ature, 9. These equations are solved with the boundary conditions shown 
in Figure 4.1. Note that the computational domain depicted in Figure 4.1 
includes an inlet entrance region that extends to 10% of the plate height H 
below the plates. This region serves to simulate in the numerical solution 
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the gradual variation of flow and temperature that would occur at the inlet 
of the actual physical problem shown in Figure 2.1. 

Recall that, for the natural convection case, the mass inflow is not known 
and has to emerge from the numerical solution. This means that the value 
of ^6 shown as a boundary condition in Figure 4.1 has to be arrived at 
through an iterative procedure. This procedure, which is described in Sec- 
tion 2.6.3 (subsection “Mass Flow Rate for Natural Convection”), needs that 
the balance of momentum, buoyancy and viscous forces, resulting in Equation 
(2.31), be satisfied. Various trial values of ^6 are used till the correct that 
yields temperature and velocity fields satisfying Equation (2.31) is obtained. 
These calculations for natural convection are thus very computationally in- 
tensive, as every solution requires a number of trial solutions to be computed. 

Although we always compute, and sometimes depict, the temperature 
and flow fields, our main concern is with the optimum spacing between the 
parallel plates which will give the highest heat transfer rate. Obviously, for 
every Grashof number, Gth, we might obtain a different optimum spacing 
or, equivalently, a different optimum aspect ratio A,.. So we need to find, 
for each Gra, the for which either the average Nusselt number Aua* 
(for the single channel case) or NuavAr (for the multiple channel case) are 
the maximum (the logic of maximizing NuavAr in the later case is given in 
Section 2.7). This requires that the Nusselt number be computed from the 
temperature field. The optimum spacing is primarily affected by the Gr^. 
We have considered a set of Grg = 10^, 10^, and 10® for our problem. 


4.1.1 Fluid Flow and Heat Transfer Patterns 

In this section, we briefly describe the velocity and temperature fields ob- 
tained in our numerical solutions. The emphasis here is not so much on the 
quantitative results but rather on the qualitative features of the solutions and 
the physics involved in obtaining them. The aim below is to demonstrate 
that the numerical solutions are quite compatible with our understanding of 
the physical processes operating in parallel plate channels. This will be done 
by considering the velocity and temperature fields for certain specific cases 
involving three different Grashof numbers, Gth = 10®, 10“^, and 10®, and four 
different aspect ratios. 
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Figure 4.1: Computational domain with boundary conditions for stream- 
function, vorticity and temperature. 
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Figure 4.2 shows the velocity vector plots for an aspect ratio, Aj. = 2.50 
and for three different Grashof numbers. In this figure, and in subsequent 
figures, the vertical and horizontal distances are shown, normalised by H, 
from the left corner of the inlet to the channel. Note that the plate starts at 
the vertical coordinate 0, the small region below being inlet upstream sec- 
tion included in the computational domain. It is evident from the figure that 
the velocity fields are roughly parabolic as expected and that the maximum 
velocity, attained at the central line, increases with Grashof number. In the 
case of the highest Grashof number, the velocity field is flatter in the central 
region indicating that the wall velocity boundary layers do not penetrate into 
the central zone of the channel. This too is expected as the boundary layer 
would be thinner for a higher Grashof number. 

Figure 4.3 shows the streamfines for the same three cases. While the 
streamlines plots themselves are not very informative, the maximum stream 
function value indicated below the plots show that the induced volume flow 
rate is a strong function of Grashof number. 

Figure 4.4 shows the isotherms for the same three cases. The difference 
between these plots is quite dramatic. The Figure 4.4(a) shows that for the 
lower Grashof number, Gr^ = 10^, the fluid entering from the bottom of the 
channel is quickly heated as it rises and attains temperature very close to 
the wall temperature even half way along the channel. Figure 4.4(b) for the 
intermediate Grashof number shows the clear presence of thermal boundary 
layers and the central line temperatures do not come close to the wall tem- 
perature even at the exit of the channel. The plot for the highest Grashof 
number, Crg — shows even a stronger influence of the thermal boundary 
layers, with the central line temperatures remaining essentially at their inlet 
values deep into the channel. 

We now consider the cases of four different aspect ratio channels for the 
intermediate Grashof number, Gr^ = 10^. The aspect ratios considered are 
1.000, 1.515, 2.000, and 3.125. It is to be noted that the higher aspect ratio 
channels are narrower. Figure 4.5 shows the velocity vector plots for these 
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Figure 4.2: Vector plots for three different Grashof numbers and aspect ratio, 
At = 2.500, for the case of natural convection: (a) Gth = 10^; (b) Grn = 
10*: (c) Grir = 10^ 
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Figure 4.4: Isotherms for three different Grashof numbers and aspect ratio, 
Ar = 2.500, for the case of natural convection: (a) Gth = 10^; (b) Gth = 
10^ (c) Grn = 10®. 
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Figure 4,5: Vector plots for four different plate-to-plate channel spacing and 
Grff = 10^, for the case of natural convection: (a) Ar = 1.000; (b) Af = 
1.515; (c) A = 2.000; (d) A, = 3.125. 
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Figure 4.6: Streamlines for four different plate- to-plate channel spacing and 
Gri{ = 10'*, for the case of natural convection: (a) Ar = 1.000; (b) Ar = 
1.515; (c) Ar = 2.000; (d) A, = 3.125. 






84 CHAPTER 4. RESULTS AND DISCUSSION 

four cases. The lowest aspect ratio channel, i.e., the widest among the fo 
shows a relatively flat profile in the central zone indicating that the widt^'^f 
the channel is much greater than the wall boundary layer thickness faith °}! 
the four figures are scaled to be equally wide, the actual channel widthr^^ 
quite different, as indicated by the horizontal coordinate). As the as 
ratio increases, the profile becomes closer to parabolic, i.e., closer to 
developed. The magnitude of the maximum velocity also decreases fortt^ 
higher aspect ratio (narrower) channels. Figure 4.6 again shows the stream^ 
lines. The volume flow rate (maximum stream function value) indicated als' 
show a decrease with increasing aspect ratio. 


Figure 4.7 shows the isotherms for the same cases. It is quite clear that 
for the lowest aspect ratio (widest) channel, the thermal boundary layers do 
not penetrate into the central zone. As the channel becomes narrower the 
boundary layers meet at the central line. For the narrowest channel [Fimre 
4.7 (d)], thermal boundary layers meet quite close to the inlet of the channel. 
Interestingly, as will be seen later, the case depicted in Figure 4.7(b) is close 
to the optimum spacing for the single channel case at this Gr^, while Figure 
4.7(d) shows the case close to the optimum spacing for the multiple channel 


These figures show the physics involved in determining the optimum spac- 
mg. It IS clear that the spacing in case (a) can be reduced further without 
significantly affecting the heat transfer as the boundary layers do not meet 
and the central zone is essentially still at the inlet temperature. However 
once the boundary layers meet, further decrease in the spacing would reduce 
the volume flow rate without much increase in the heat transfer per unit 
volume flow as the outflow temperature is already close to the maximum 
attmnable So the optimum spacing for the single channel case seems to be 
that which allows the thermal boundary layer to meet somewhere close to 
the outflow plane. For the multiple channel case, we maximize and 

can thus decrease the spacing further as the resulting increase in the aspect 
ratio compensates for the decrease in up to a point. Thus, optimum 

pacing for the multiple channel case occurs at a higher aspect ratio as shown 
in figure 4.7(d). 
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Figure 4.7: Isotherms for four different plate-to-plate channel spacing and 
Gr[j = 10'*, for the case of natural convection: (a) Ar = 1.000; (b) Ar = 
1.515; (c) Ar = 2.000; (d) = 3.125. 
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Table 4.1 Dimensionless induced flow rate(A^), average Nussplt . . 

and the product of average Nusselt number and aspect ratin 

aspect ratio of the channel at Pr = 0.70, Grfj = 10^ ^ ‘different 


1 ^ 


NxLq^v 

NUa. 

•^Uav-^r 

av 

1 0.556 

1.293 

3.688 

3.695 

2.049 

2^0^3^ 

0.629 

1.126 

3.707 

3.691 

2.331 

2.321 1 

0.667 

1.054 

3.715 

3.689 

2.476 

2.459 

0.680 

1.030 

3.718 

3.689 

2.529 

2.509 1 

0.694 

1.006 

3.721 

3.688 

2.584 

2.561 

0.741 

0.934 

3.729 

3.687 

2.763 

2.731 

0.775 

0.886 

3.735 

3.686 

2.895 

2.857 

0.813 

0.838 

3.740 

3.684 

3.041 

2.995 

1 0.833 

0.814 

3.743 

3.683 

3.119 

3.069 

0.855 

0.792 

1 3.748 

3.745 

' 3.203 

' 3.201 

1.000 

0.651 

3.677 

3.650 

3.677 

3.650 

1.333 

0.449 

3.631 

3.571 

4.841 

4.762 

1.786 

0.285 

3.055 

3.002 

5.455 

5.360 

1.852 

0.267 

2.952 

2.892 

5.466 

5.355 

1.923 

0.249 

2.838 

2.771 

5.458 

5.329 

2.000 

0.231 

2.697 

2.638 

5.393 

5.276 

2.250 

2.500 

0.181 

0.141 

2.298 

1.953 

2.217 

1.897 

5.171 

4.882 

4.989 

4.743 

3.333 

0.066 

1.320 

1.358 

4.400 

4.526 


4.1.2 Optimum Spacing and Voiume Piow Rate in the 
Natural Convection Case 

Tabte rf 4 2 “‘"“I “ 

each caae varies somewhat and is cht^en so t “ 

be obtained for the chosen Gr„. ^ optimum spacing can 


The average Nusselt number, Nu 
has been computed in two different 


av over the heated walls of the channel 
ways so as to provide a cross-check of 
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Table 4.2 Dimensionless induced flow rate (AT^), average Nusselt number, 
and the product of average Nusselt number and aspect ratio for different 
aspect ratio of the channel at Pr = 0.70, Gth = 10'^ 


Ar 


N Uav 



NUavAr 

1.000 

0.747 

6.278 

6.349 

6.278 

6.349 

1.250 

0.580 

6.300 

6.368 

7.875 

7.960 

1.429 

0.495 

6.353 

6.395 

9.076 

9.101 

1.471 

0.478 

6.354 

6.394 

9.345 

9.402 

1.515 

0.461 

6.355 

6.392 

9.629 

9.684 

1.563 

0.444 

6.355 

6.390 

9.930 

9.984 

1.667 

0.410 

6.352 

6.381 

10.586 

10.635 

2.000 

0.324 

6.291 

6.306 

12.581 

12.611 

2.250 

0.276 

6.174 

6.179 

13.892 

13.904 

2.500 

0.237 

5.962 

5.977 

14.906 

14.943 

2.941 

0.184 

5.411 

5.432 

15.914 

15.977 

3.125 

0.166 

5.125 

5.155 

16.017 

16.110 

3.333 

0.148 

4.791 

4.816 

15.968 

16.054 

3.571 

0.133 

4.318 

4.340 

15.422 

15.499 

4.000 

0.105 

3.561 

3.572 

14.245 

14.287 

5.000 

0.061 

2.183 

2.178 

10.914 

10.890 
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Table 4.3 Dimensionless induced flow rate (AT^) averaee Wucc i+ 
and the product of average Nusselt mimKoT- o J Nusselt number, 

aspect ratio of the channel at Pr = 0.70, Gvff = different 


Ar 


Nudv 

1.333~ 

0.597 

11.343 

1.389 

0.570 

11.342 

1.429 

0.552 

11.343 

1.449 

0.544 

11.343 

1.471 

0.543 

11.335 

1.515 

0.521 

11.337 

1.613 

0.488 

11.339 

1.667 

0.468 

11.311 

2.000 

0.378 

11.298 

2.500 

0.287 

11.268 

3.333 

0.197 

11.204 : 

5.000 

0.114 

9.918 

5.263 

0.104 

9.531 

5.555 

0.095 

9.071 

5.882 

0.086 

8.513 

6.667 

0.067 1 

7.086 1 


Nu„ 

av 

11.214 

11.213 

11.215 

11.214 
11.205 
11.208 
11.210 
11.189 
11.218 
11.196 
11.109 
9.855 
9.439 
8.974 
8.412 
6.989 


W A 
15.123 
15.753 
16.204 
16.439 
16.669 
17.178 
18.288 
18.852 
22.595 
28.171 
37.347 
49.588 
50.162 
50.393 
50.076 
47.237 


14.952 

15.574 

16.021 

16.252 

16.478 

16.982 

18.080 

18.649 

22.436 

27.990 

37.030 

49.275 

49.680 

49.852 

49.483 

46.594 
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the accuracy of the numerical results. The one which is based on the normal 
temperature gradient at the wall is 



The second method of computing the average Nusselt number is based on the 
energy balance of the heat flux from the walls and at the inlet and the outlet 
of the channel. For the latter, both the convective and conductive fluxes are 
considered. This balance yields the following dimensionless equation for the 
average Nusselt number: 


=l{_-PT^sJ'U6\x,,dY-jA§lx,idY (4.2) 

+ /o* §i lx=o dY + fA Ue\x,i dY} 

where the first term on the right corresponds to convective heat flux at the 
outlet, the second term to conductive heat flux at the inlet, the third term 
to conductive heat flux at the outlet, and the fourth term to convective heat 
flux at the inlet of the channel. 

Special care was taken for computing the second and the fourth integral 
of These integrals were actually computed along a surface slightly 

outside the channel entrance in order to avoid the expected leading-edge sin- 
gularity at the inlet of the channel. It was noticed that while the conduction 
term at the outlet of the channel was negligibly small (as the flow there is 
convection dominated), the outward (downward) flow of heat at the inlet due 
to conduction is surprisingly significant . 

The two alternative ways of computing the Nusselt number provide us 
with a means of cross-checking the accuracy of the numerical simulations. 
It can be seen from Tables 4.1, 4.2, and 4.3 that the diflference between the 
two values of the Nusselt numbers is never more than 3.5% and usually lies 
within 1.5%. This degree of correspondence is reassuring. In the plots that 
follows, only Nuav (Equation 4.1) is presented, as the two different Nusselt 
numbers would be indistinguishable on a graph. 
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Table 4.4: Comparison of the present work with Anand et al. (1992), with 
rega rds to maximum average Nusselt number 


Grashof number, 
Grn 

Maximum average Nusselt number iNuav)mn-r 

Anand et al, 1992 

Present study 

10® 

3.55 

3.74 

10^ 

6.15 

6.36 

10® 

10.85 

11.34 


Single Channel Case 

The Figure 4.8 shows the average Nusselt number as a function of the as- 
pect ratio for all three Grashof numbers. It can be seen that, as expected 
the Grashof number has a significant effect on the Nuav and hence the heat 
transfer. For the single channel case, the maximum Nua^ would determine 
the optimum spacing. It can be seen from the figure that such an optimum 
does exist. This is because as the aspect ratio increases (i.e., the width de- 
creases), the flow rate decreases but the temperature gradients increase and 
the overall heat transfer may be increased. As the aspect ratio increases, the 
“chimney effect” would also create a greater draft and hence greater heat 
transfer. These contrasting effects ensure that an optimum plate-to-plate 
spacing does exist which is greater than zero. However, from the Figure 4.8 
we can see that the N u^v dependency on aspect ratio is minimal for low Ar. 
So the optimum is not a strong one. In fact, it can be argued that any aspect 
ratio below 0.85 for Gth = 10^ 1.67 for Gth = 10^ and 2.00 for Gru = 10® 
would be close to optimum for the single chaimel case. 

Table 4.4 shows the comparison of the present results with that of Anand 
et al. (1992) who solved the same problem using the boundary layer equa- 
tions. The table shows the maximum Nuan obtained in the two studies. It 
can be seen that the two results are quite close, thereby showing that the 

boundary layer approximation is surprisingly accurate for the quantities con- 
sidered. 
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Table 4.5 Range of optimum aspect ratios and maximnm ^ 

for the single and multiple channel cases ^ ^ transfer rate 


Grff 

Single chi 
i^r)opt 

innel case 

^^'^av)max 

Multiple 

lU 

1 

0.81 - 0.85 

3.74 

1.79 - 1.92 

5 46 

i(J 

1 a 5 

1.25 - 1."^ 

n 6.36 

2.94 - 3.33 

^ 1®~ 

iU 

1.33 - 2.00 

11.34 

5.26 - 5.88 

5S~ 


Multiple Channel Case 

For the multiple channel case, as was argued in Section 9 7 i-h 

of the product Nn„A^ determtoes the optimum A Fi^re fo T'T 

quantity vemus aspect ratio, 4 for the ftree^li c 

(Grn = ^0^ 10<, and 10»). It may be seen that the optSum k “T ® 

demarcated in fee cases as compared to those ina^^nSe 

Table 4.5 shows the optimum spacing and the values Nu A T' 

qmte evident. As Gr„ increases, optimum aspect rltfc “ 

the optimum channel width decreases. ^ ™ increases, i.e., 

As will be seen below, when Crtr is kpnt mTicfc.Ta+ ^ ^ a ■ . 
volume fiow rate drops quite sharply leading to lower heaUw'"^^' 
Gru IS mcteased, higher volume flow rate and hence heat 
be mamfeed for a greater range of aspect ratios. Th^fom 
4 tends to mcrease with increasing Gr„ Practicallv °P*™™ 

Ugber Gr„ the channel walls shoufd o^ima^tellosL “““ 


Volume Flow Rate 

Fi^e 4.10 shows the induced flowrate, A4>, values for different 4 n 
Thm quantity gives the nondimensional volume flow rfe A® f"' 

ivhere, u is the volnme flow rate per unit depth of the diiineT 7t 

seen m Firare 4 10 that fhn . i i ^ cnannel. It can be 

slightly wiScrX a co-^^^^^ 

6 biy wiin LtTh tor a constant aspect ratio. However this is nniv or. r ; 
of the nouamensionalisation used. Figm. 4.11 sho« the ssTe dl 
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Figure 4.9: Variation of the product of average Nusselt number and aspect 
ratio with aspect ratio at three different Grashof numbers, for the natural 
convection case. 
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Figure 4.10: Variation of din 
ratio at three different Gras) 
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different nondimensionalisation of the volume flow rate, A^* = ^. It can be 
seen that the effect of increasing Grn is quite dramatic, the volume flow rate 
increasing dramatically with Grn- As the computations were not done for 
very low aspect ratios (i.e., large channel width), the “chimney effect” is not 
discernible in these figures, for the flow rate decreases monotonically with 
increasing aspect ratio, in all cases. 


4.2 Mixed Convection Case 

In this section, the results for the mixed convection case are presented for 
both the single and the multiple channel cases for Pr = 0.70, Gth = 10®, 10^, 
and 10®, and buoyancy parameter values of =0.1, 0.5, 1.0, and 1.5. 

The physical situation is the same as depicted in Figure 2.1, except that 
an external flow is superimposed on the buoyancy-induced natural convec- 
tion. The governing equations solved for this case are (2.26), (2.29), and 
(2.30) which use the false transient method to obtain the steady state so- 
lutions. These equations are solved with the boundary conditions shown in 
Figure 4.1. Note that the stream function, vorticity and the velocity compo- 
nents are nondimensionalised differently than in the natural convection case. 
Furthermore, unlike in the natural convection case, the stream function value 
at the right plate, is known a priori and is given by 'Fj = —Ar~^ where, 
Ar is the aspect ratio of the channel (the dimensional ■06 = —UcL, nondi- 
mensionalised by dividing by UcH, yields the non-dimensional is -L/H 
or — Ar"^). With this and the other boundary conditions, the problem can 
be solved by the usual method without iterating for the volume flow rate as 
was done in the natural convection case. 

4.2.1 Fluid Flow and Heat Transfer Patterns 

In this section, we qualitatively describe the numerical solutions to show 
that they are compatible with the physics of flows in vertical parallel plate 
channels. We consider specific cases involving three different Grashof num- 
bers, four different aspect ratios, and four different values of the buoyancy 
parameter. 

Figure 4.12 shows the velocity vector plots for an aspect ratio, At = 2.50 
and for three Grashof numbers, Gr/r = 10®, 10^, and 10® for a fixed buoyancy 
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parameter of 1.0. It is apparent from the figure that the velocity profiles are 
roughly parabolic for the lower Grashof numbers indicating that the flow 
becomes fully developed quickly along the channel. However, for the high- 
est Grashof number, the velocity profile in the central region of the channel 
is relatively flat indicating that the boundary layers have not reached the 
central core of the channel. Figure 4.13 shows the streamlines for the same 
cases. The corresponding non-dimensional volume flow rate, (=vfUcH 
= where v is the volume flow rate per unit depth; Recall that [see 

Chapter 2] the dimensionless variables are differently defined for the natural 
and mixed convection cases) is the same in all cases as Ar is kept constant. 
However, the dimensional volume flow rate increases with Grashof number, 
for a fixed buoyancy parameter. However, unlike in the natural convection 
case, the inlet velocity in this case is essentially imposed and the volume flow 
rate increases only because a higher Grashof number for a fixed buoyancy 
parameter implies a higher Reynolds number, and thus a greater flow rate. 

Figure 4. 14 shows the isotherms for the same three cases. The differences 
between these plots are quite apparent. For the highest Grashof number, Gru 
(and Rch) case, the flow is obviously not fully developed and the boundary 
layers are clearly evident. For the low Gth (and Reu) cases the bound- 
ary layers meet in the flow domain and the outlet temperature is quite close 
to the plate temperature, particularity so for the lowest Gth (and Rea) case. 

Figure 4.15 shows the velocity vector plots for Gvh = 10^ and = 1.0, 
and four different aspect ratios 1.6667, 2.7778, 3.3333, and 5.0000. Here, 
the flow tends to become fully developed for the highest aspect ratio (i.e., 
the narrowest channel) while remaining relatively flat near the central line 
for the lowest aspect ratio case. Figure 4.16 shows the streamlines for the 
same cases indicating the corresponding non-dimensionalised volume flow 
rate, A^. Here, as both the Grashof munber and buoyancy parameter are 
fixed, the volume flow rate (both dimensional and non-dimensional) is seen 
to decrease with increasing aspect ratio, as expected. 

Figure 4.17 shows the isotherms for the same cases. It can be seen that 
the thermal boundary layers tend to meet closer to the inlet as the aspect 
ratio is increased. Interestingly, the case shown in Figure 4. 17(b) corresponds 
to the optimum spacing (i.e. with the highest Nuav) for the single channel 
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Figure 4 J2^Vector plots for three different Grashof numbers and Ar = 2.50, 

convection: (a) Crg = 10^; (b) 
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Figure 4.13: Streamlines for three different Grashof numbers and Ar = 2.50, 
Crn/Re^j = 1.0, for the case of mixed convection: (a) Grn = 10^; (b) 
Grij = 10^ (c) Grn = 10^ 
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Figure 4.14: Isotherms for three different Grashof numbers and At = 2.50, 
GrH/R£% — 1-0) for the case of mixed convection: (a) Grn = (^) 

Gth = lO"*; (c) Gtj{ = 10^. The volume flow rate =0.4 for all cases. 
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Figure 4.15: Vector plots for four different channel spacing and Grn = 10^, 
CrHlRCfi = 1-0, for the case of mixed convection: (a) Ar = 1.6667; (b) 
Ar = 2.7778; (c) Ar = 3.3333; (d) Ar = 5.0000. 
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t ? = ‘O*' 

n fi- rh) A convection: (a) Ar = 1.6667, A$ = 
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(c) (d) 


Figure 4.17: Isotherms for four different channel spacing and Crg = 10'^, 
GTjjjRe\ = 1.0, for the case of mixed convection: (a) A^. = 1.6667; (b) 
.4, = 2.7778; (c) A, = 3.3333; (d) A^ = 5.0000. 
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at this Gth = 10^. As the aspect ratio is increased the channel become 
narrower and the thermal gradient should become steeper. However, the 
volume flow rate of the imposed flow also decreases. These counter-acting 
effects seem to yield an optimum Nuav when the boundary layers meet close 
to half way down the channel. For the multiple channel case, the net flow 
rate per unit width does not decrease as aspect ratio increases. Hence, the 
product NUavAr keep on increasing with aspect ratio finally to an asymp- 
totic maximum value at high aspect ratio (this will be shown later). Figure 
4.17(d) shows a case very close to this asymptotic maximum. 

Figure 4.18 shows the velocity vector plots for four different buoyancy pa- 
rameters (= 0.1, 0.5, 1.0, and 1.5) at fixed Ar = 2.5 and Grn = 10^. It can 
be seen that as the buoyancy parameter increases from 0.1 to 1.5 (and Ren 
decreases) the velocity boundary layers become thicker as expected. While 
the velocity profile at = 0.1 is flat throughout the central core of the 
diannel, it tends to become parabolic as buoyancy parameter increases (and 
Ren decreases). Figure 4.19 depicts the streamlines for the same conditions 
showing corresponding dimensionless volume flow rate. Although the nondi- 
mensional volume flow rate remains constant for all the cases shown, as the 
buoyancy parameter increases and Ren decreases, the dimensional volume 
flow rate decreases. 

Figure 4.20 illustrates isotherms for four different buoyancy parameters 
( = 0.1, 0.5, 1.0, and 1.5) and Gth = 10^. It is clearly evident from these 
isotherms that as the buoyancy parameter increases (Res decreases), the 
thermal boundary layer becomes thicker. At = 0.1 (i.e., the highest 

Ren), the two thermal boundary layers do not meet inside the channel, indi- 
cating that inter-plate spacing has to be decreased further in order to obtain 
an optimum spacing. As ^ increases to 0.5 and above, the two thermal 

boundary layers meet inside the channel and this meeting point moves to- 
wards the inlet of the channel. 


4.2.2 Optimum Spacing in the Mixed Convection Case 

The numerical results are tabulated for the mixed convection case in Tables 
4.6-4.17 for Gth = 10^ 10^, and 10® and gf = 0.1, 0.5, 1.0, and 1.5. The 

H 
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Figure 4.18: Vector plots for four different channel buoyancy parameters and 
Grff = 10^, Ar = 2.50, for the case of mixed convection: (a) GrnlRe^H — 0-1; 
(b) GrnfRel = 0.5; (c) GthIRc^ = 1.0; (d) GralRe^ = 1.5. 
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font' 2 loinT ‘T ^ 

OrWRel = oi; W W =To ^ 

rate = o.40 for all cases. ^ - 1.5. The v< 


parameters and Grjf = 
0 GrH/Re% = 0.1; (b) 
= 1.5. The volume flow 
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(c) (d) 


Figure 4.20: Isotherms for four different buoyancy parameters and Gru = 
10^, Ar = 2.50, for the case of mixed convection: (a) GrulRe^ = 0.1; (b) 
Orn/Rel = 0.5; (c) Gth/RcI = 1.0; (d) GralRel = 1.5. 
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Table 4.6 Variation of pressure drop, average Nusselt number and 
product of average Nusselt number and aspect ratio, with aspect raK r 
mixed convection case at Pr = 0.70, Grn = 10® and ^ = n i 


1 A 

APa« 


NUa. 



1.000 

0.276 

6.745 

6.702 

6.745 

6.702 

2.000 

0.322 

7.374 

7.382 

14.747 

14.763 

2.500 

0.353 

7.600 

7.621 

18.999 

19.051 

2.778 

0.374 

7.667 

7.679 

21.296 

21.330 

3.030 

0.396 

7.677 

7.681 

23.262 

23.276 

3.125 

0.404 

7.660 

7.669 

23.937 

23.964 

3.226 

0.413 

7.677 

7.650 

24.764 

24.677 

3.333 

0.423 

7.650 

7.622 

25.500 

25.406 

3.448 

0.434 

7.613 

7.584 

26.253 

' 26.150 

3.571 

0.446 

7.559 

7.532 

26.997 

26.902 

5.000 

0.592 

6.413 

6.378 

31.889 

32.064 

6.667 

0.786 

5.037 

5.058 

33.582 

33.719 

7.143 

0.842 

4.798 

4.739 

34.272 

33.851 

7.692 

0.908 

4.392 

4.412 

33.782 

33.937 

8.333 

0.985 

4.057 

4.078 

33.811 

33.981 

10.000 

1.185 

3.380 

3.400 

33.804 

.. 

34.004 


Pmdtl nmnber, m ^ these cases is 0.70. The average Nusselt numbe.s 
have been computed as before in two different ways; Nu., using the tem- 
Pfflaluiedenvative at the wall and JVu„ by using the heat balLe. The 
diff« bet^ these two computed values are generally within 2 to 3 per 
cent. Thus only JVuqp is presented in the figures that follow. 


Single Channel Case 

In the single channel ca^ the optimum aspect ratio is determined simply 
by the m^mum ATu,.. Figures 4.21, 4.22, and 4.23 show the average Nus- 

Or “"l«f To-'Z "“"“‘TV “<1 

10 , , and 10 respectively. It can be seen that, unlilie in the nat- 

ural convection case, the maxima here are quite clear. However, it can also 
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Table 4.7 Variation of pressure drop, average Nusselt number, and the 
product of average Nusselt number and aspect ratio, with aspect ratio for 
mixed convection case at Pr = 0.70, Gru = 10^, and = 0.5 

H 


Ar 

APa„ 

NUav 




1.667 

0.444 

5.227 

5.298 

8.712 

8.831 

1.754 

0.459 

5.246 

5.312 

9.204 

9.319 

1.786 

0.470 

5.220 

5.292 

9.322 

9.450 

1.852 

0.477 

5.257 

5.317 

9.736 

9.846 

1.923 

0.495 

5.229 

5.293 

10.056 

10.178 

1.961 

0.497 

5.258 

5.311 

10.310 

10.414 

2.000 

0.509 

5.225 

5.285 

10.451 

10.570 

2.500 

0.613 

5.039 

5.098 

12.599 

12.744 

3.333 

0.815 

4.441 

4.462 

14.804 

14.873 

5.000 

1.253 

3.093 

3.127 

15.463 

15.636 

5.263 

1.325 

2.946 

2.977 

15.507 

15.668 

5.556 

1.406 

2.798 

2.825 

15.546 

15.692 

5.882 

1.496 

2.649 

2.671 

15.585 

15.710 

6.250 

1.597 

2.499 

2.514 

15.619 

15.713 

6.667 

1.711 . 

2.291 

2.359 

15.275 

15.725 

7.143 

1.842 

2.137 

2.201 

15.266 

15.724 
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Table 4.8; Variation of pressure drop, average Nusselt number, and the prod- 
uct of average Nusselt number and aspect ratio, with aspect ratio for mixed 
convection case at Pr = 0.70, Gvfj = 10^, and = 1.0 


Ar 

APav 

■A^ '^av 

Nu,, 

'lidyAj' 

NUavAr 

1.449 

0.413 

4.511 

4.601 

6.538 

6.668 

1.515 

0.430 

4.518 

4.602 

6.845 

6.973 

1.587 

0.449 

4.519 

4.598 

7.173 

7.299 

1.667 

0.470 

4.513 

4.587 

7.521 

7.645 

1.754 

0.495 

4.498 

4.566 

7.891 

8.011 

1.852 

0.524 

4.471 

4.534 

8.280 

8.396 

1.961 

0.557 

4.430 

4.486 

8.686 

8.796 

2.000 

0.572 

4.364 

4.453 

8.729 

8.907 

2.083 

0.596 

4.368 

4.417 

9.099 

9.203 

2.222 

0.644 

4.241 

4.319 

9.424 

9.597 

2.500 

0.740 

4.044 

4.110 

10.110 

10.274 

4.000 

1.309 

2.791 

2.860 

11.163 

11.438 

5.000 

1.711 

2.224 

2.318 

11.120 

11.588 

5.882 

2.064 

1.898 

1.981 

11.162 

11.652 

6.667 

2.376 

1.679 

1.753 

11.192 

11.684 

7.143 

2.564 

1.568 

1.637 

11.199 

11.694 
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Table 4.9: Variation of pressure drop, average Nusselt number, and the prod- 
uct of average Nusselt number and aspect ratio, with aspect ratio for mixed 
convection case at Pr = 0.70, Gth = 10^, and = 1.5 


Ar 

XPav 

^ '^av 


^ '^av-^T 

NUaAr 

1.333 

0.322 

4.161 

4.261 

5.547 

5.681 

1.389 

0.339 

4.162 

4.257 

5.780 

5.913 

1.449 

0.359 

4.158 

4.249 

6.027 

6.158 

1.515 

0.382 

4.150 

4.236 

6.288 

6.419 

1.587 

0.407 

4.135 

4.216 

6.563 

6.692 

1.667 

0.436 

4.111 

4.187 

6.852 

6.978 

1.754 

0.469 

4.077 

4.148 

7.153 

7.277 

1.852 

0.507 

4.031 

4.095 

7.464 

7.584 

1.961 

0.551 

3.969 

4.027 

7.782 

7.896 

2.083 

0.603 

3.887 

3.938 

8.098 

8.204 

2.500 

0.793 

3.524 

3.592 

8.809 

8.981 

3.333 

1.215 

2.769 

2.863 

9.230 

9.542 

5.000 

2.039 

1.876 

1.974 

9.379 

9.870 

5.882 

2.481 

1.606 

1.692 

9.446 ' 

9.952 

6.667 

2.868 

1.425 

1.500 

9.497 

9.999 

7.143 

3.101 

1.332 

1.402 

9.514 

10.016 

7.692 

3.368 

1.242 

1.307 

9.552 

10.053 
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Table 4.10 Variation of pressure drop, average Nusselt number, and the 
product of average Nusselt number and aspect ratio, with aspect ratio for 
mixed convection case at Pr = 0.70, Grn = 10'^, and = 0.1 


Ar 

APa„ 

N Uav 

Nul 

•AT lifiyAf 

mmm 

1.667 

0.145 

11.585 

11.468 

19.308 

■feUhai 

1.961 

0.149 

11.832 

11.729 

23.199 

22.998 

2.222 

0.155 

12.327 

12.046 

27.394 

26.769 

2.500 

0.158 

12.285 

12.126 

30.712 

KtMMII 

3.333 

0.171 

13.071 

12.767 

43.570 

42.557 

5.000 

0.200 

13.358 

13.308 

66.789 

66.540 

5.263 

0.208 

13.374 

13.314 

70.388 

70.073 

5.556 

0.217 

13.359 

13.288 

74.217 

73.824 

5.882 

0.227 

13.304 

1 

13.221 

78.257 

77.771 

6.667 

0.253 

13.017 

12.906 

86.777 

86.041 

7.692 

0.288 

12.381 

12.253 

95.239 

94.253 

10.000 

0.372 

10.513 

10.382 

105.125 

103.824 

12.500 

0.466 

8.637 

8.526 

107.959 

106.580 

16.667 

0.624 

6.502 

6.424 

108.368 

107.071 

1 20.000 

0.751 

5.415 

5.355 

108.307 

107.095 
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Table 4.11 Variation of pressure drop, average Nusselt number, and product 
of average Nusselt number and aspect ratio with aspect ratio for mixed 
convection case at Pr = 0.70, Gth = 10^, and ^ = 0.5 

H 


Ar 

^Pav 

iV U(xy 

NUav 


NUavAr 

1.667 

0.178 

8.368 

8.381 

13.947 

13.969 

1.961 

0.191 

8.581 

8.586 

16.826 

16.835 

2.222 

0.207 

8.847 

8.783 

19.659 

19.519 

2.500 

0.219 

8.984 

8.935 

22.459 

22.337 

2.857 

0.237 

9.041 

9.019 

25.831 

25.768 

3.125 

0.252 

9.109 

9.094 

28.465 

28.418 

3.333 

0.264 

9.202 

9.125 

30.673 

30.415 

3.571 

0.278 

9.192 

9.124 

32.829 

32.587 

3.704 

0.287 

9.178 

9.110 

33.992 

33.740 

4.000 

0.306 

9.110 

9.041 

36.439 

36.163 

5.000 

0.373 

8.318 

8.320 

41.588 1 

41.602 

6.667 

0.510 

6.985 

6.969 

46.568 

46.458 

10.000 

0.800 

4.805 

4.795 

48.050 

47.947 

11.111 

0.897 

4.322 

4.318 

48.023 

47.982 

12.500 

1.018 

3.843 

3.837 

48.038 

47.967 
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Table 4.12; Variation of pressure drop, average Nusselt number, and the 
product of average Nusselt number and aspect ratio, with aspect ratio for 
mixed convection case at Pr = 0.70, Gth = 10^, and = 1.0 

M 


Ar 

lAPav 

^ '^av 

NUav 

^ ^aV’^r 

Nu^^Ar 

1.754 

0.157 

7.404 

7.435 

12.989 

13.044 

1.961 

0.172 

7.520 

7.542 

14.745 

14.788 

2.222 

0.196 

7.692 

7.670 

17.092 

17.045 

2.500 

0.216 

7.764 

7.748 

19.411 

19.370 

2.564 

0.221 

7.777 

7.760 

19.941 

19.896 

2.778 

0.239 

7.800 

7.779 

21.666 

21.608 

2.857 

0.245 

7.726 

7.733 

22.075 

22.094 

3.030 

0.261 

7.788 

7.763 

23.599 

23.524 

3.333 

0.289 

7.681 

7.648 

25.602 

25.494 

5.000 

0.465 

6.353 

6.380 

31.764 

31.898 

5.556 

0.533 

5.892 

5.914 

32.732 

32.855 

6.667 

0.674 

5.038 

5.057 

33.586 

33.714 

8.333 

0.889 

4.059 

4.077 

33.826 ! 

33.976 

9.091 

0.986 

3.722 

3.740 

33.840 

33.999 

10.000 

1.102 

3.382 

3.400 

33.823 

33.002 
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Table 4.13 Variation of pressure drop, average Nusselt number, and the prod- 
uct of average Nusselt number and aspect ratio, with aspect ratio for mixed 
convection case at Pr = 0.70, Grn = 10'^, and = 1.5 

H 


Ar 

APa. 

NU(iy 




1.961 

0.125 

6.955 

6.987 

13.638 

13.700 

2.222 

0.154 

7.067 

7.058 

15.705 

15.685 

2.381 

0.170 

7.079 

7.075 

16.855 

16.845 

2.564 

0.188 

7.083 

7.076 

18.161 

18.144 

2.778 

0.212 

7.057 

7.046 

19.601 

19.573 

3.030 

0.241 

6.985 

6.970 

21.167 

21.123 

3.333 

0.278 

6.847 

6.827 

22.824 

22.757 

5.000 

0.512 

5.332 

5.374 

26.658 

26.870 

6.667 

0.783 

4.128 

4.164 

27.520 

27.760 

7.143 

0.861 

3.858 

3.892 

27.560 

27.802 

7.692 

0.951 

3.544 

3.618 

27.569 

27.828 

8.333 

1.055 

3.310 

3.341 

27.581 

27.844 

9.091 

1.177 

3.033 

3.064 

27.577 

27.851 

10.000 

1.321 

2.756 

2.785 

27.562 

27.849 
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Table 4.14 Variation of pressure drop, average Nusselt number, and the prod- 
uct of average Nusselt number and aspect ratio, with aspect ratio for mixed 
convection case at Pr = 0.70, Gth = 10®, and^f = 0.1 

H 


Ar 


NU(iy 


^ '^av^r 


2.222 

0.076 

21.136 

19.752 

46.970 

43.893 

3.333 

0.082 

21.464 

20.648 

71.545 

68.827 

4.000 

0.083 

22.119 

21.013 

88.475 

84.051 

5.000 

0.085 

22.259 

21.785 

111.297 

108.924 

6.667 

0.091 

23.239 

22.843 

154.929 

152.284 

7.692 

0.100 

23.661 

23.282 

182.005 

179.091 

8.333 

0.106 

23.858 

23.468 

198.817 

195.563 

9.091 

0.113 

23.974 

23.570 

217.949 

214.270 

10.000 

0.122 

23.940 

23.521 

239.400 

235.214 

11.111 

0.134 

23.653 

23.224 : 

262.813 

258.042 

12.500 

0.149 

22.971 

22.550 

287.139 

281.877 

16.667 

0.196 

19.752 

19.374 

329.198 

322.900 

25.000 

0.295 

13.741 

13.496 

343.530 

337.411 

33.333 

0.395 

10.269 

10.119 

342.287 

337.287 
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Table 4.15 Variation of pressure drop, average Nusselt number, and the prod- 
uct of average Nusselt number and aspect ratio, with aspect ratio for mixed 
convection case at Pr = 0.70, Gr^ = 10®, and = 0.5 


A 

APa. 

NUav 

NUa. 

NUav^r 


2.381 

0.090 

14.738 

14.328 

35.091 

34.115 

3.333 

0.104 

15.311 

15.026 

51.035 

50.086 

3.704 

0.110 

15.540 

15.283 

57.556 

56.603 

5.000 

0.119 

15.761 

15.380 

78.803 

76.902 

5.882 

0.137 

15.969 

15.856 

93.938 

93.267 

6.250 

0.146 

15.997 

15.855 

99.978 

99.091 

6.667 

0.155 

15.961 

15.808 

106.410 

105.384 

7.692 

0.179 

15.645 

15.465 

120.342 

118.962 

10.000 

0.237 

14.114 

13.913 

141.143 

139.128 

12.500 

0.304 

12.040 

11.855 

150.500 

148.189 

16.667 

0.419 

9.213 

9.074 ^ 

153.542 

151.226 

20.000 

0.511 

7.684 

7.573 

153.687 

151.453 

25.000 

0.649 

6.140 

6.061 

153.499 

151.512 
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Table 4.16 Variation of pressure drop, average Nusselt number, and the prod- 
uct of average Nusselt number and aspect ratio, with aspect ratio for mixed 
convection case at Pr = 0.70, Gth = 10^, and = 1.0 

H 


Ar 

APa. 

NUfiy 

mam 


mmm 

2.222 

0.066 

12.849 

12.536 

28.552 

27.858 

3.030 

0.085 

13.207 

13.002 

40.022 

39.400 

3.333 

0.092 

13.363 

13.174 

44.542 

43.914 

4.546 

0.110 

13.397 

13.376 

60.896 

60.799 

4.762 

0.116 

13.432 

13.402 

63.960 

63.818 

5.000 

0.123 

13.452 

13.413 

67.259 

67.065 

5.263 

0.132 

13.451 

13.401 

70.795 

70.531 

5.556 

0.141 

13.421 

13.359 

74.561 

74.217 

5.882 

0.151 

13.352 

13.277 

78.539 

78.102 

6.667 

0.178 

13.042 

12.937 

86.944 ^ 

86.248 

10.000 

0.304 

10.513 

10.383 

105.130 

103.835 

12.500 

0.405 

8.636 

8.526 

107.955 

106.579 

16.667 

0.575 

6.503 

6.424 

108.385 

107.072 

1 20.000 

0.708 

5.417 

5.355 

108.336 

107.096 
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Table 4.17: Variation of pressure drop, average Nusselt number, and the 

product of average Nusselt number and aspect ratio, with aspect ratio for 

mixed convection case at Pr = 0.70, Gth = 10^, and^^ = 1.5 




Ar 

XPav 

^ '^av 

NUa. 

^ ^CLV Ar 

Nu^^Ar 

2.778 

0.046 

12.120 

11.943 

33.666 

33.174 

3.030 

0.056 

12.224 

12.070 

37.042 

36.576 

3.333 

0.065 

12.345 

12.197 

41.152 

40.655 

4.167 

0.077 

12.162 

12.180 

50.674 

50.749 

4.348 

0.084 

12.178 

12.188 

52.949 

52.992 

4.546 

0.091 

12.184 

12.185 

55.380 

55.385 

4.762 

0.099 

12.173 

12.164 

57.968 

57.926 

5.000 

0.109 

12.142 

12.120 

60.711 

60.602 

6.667 

0.180 

11.383 

11.303 

75.890 

75.353 

8.333 

0.261 

10.067 

9.974 

83.894 

83.113 

10.000 

0.344 

8.712 

8.620 

87.115 

86.197 

12.500 

0.474 

7.059 

6.987 

88.232 

87.343 

14.286 

0.566 

6.185 

6.120 

88.352 

87.430 

16.667 

0.687 

5.300 

5.247 

88.313 

87.446 
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Table 4.18 Range of optimum aspect ratio and corresponding average 
Nusselt number for different Grashof numbers and buoyancy parameters- 
Single Channel case 



^ = 0.1 

, 


0.5 


1.0 


iT^ 

Grn 

Ar 


Ar 

^ '^dv 

Ar 


4 


10^ 

2.78-3.23 

7.66 

1.75-1.96 

5.25 

1.45-1.67 

4.52 

1.33-1.45 

4.16“ 

10^ 

5.00-5.55 

13.37 

3.13-3.70 

9.20 

2.50-3.03 

7.80" 

2.22-2.78 


10® 

8.33-10.00 

23.97 

5.88-6.67 

15.99 

4“.76-5.56 

13.45 

3.03-4.55 

12.35“ 


be seen that a range of aspect ratios give values of close enough to the 
maximum that all may be considered optimum. 

Table 4.18 shows the range of optimum aspect ratio and corresponding 
Nuav for each Grn and for all buoyancy parameters considered. Some trends 
are immediately evident: as Crg increases, for a fixed buoyancy parameter, 
the optimum Ar increases (implying smaller plate-to-plate spacing) and the 
maximum Nuav also increases. While it is clear that the Nuav should in- 
crease with increasing Grir, it is not so obvious why the optimum spacing 
should decrease. However, it should be remembered that a fixed buoyancy 
parameter and increased Gth means a higher Reynolds number, and hence 
thinner momentum boundary layers. As Pr is fixed, the thinner momentum 
boundary layers directly imply thinner thermal boundary layers. As the op- 
timum spacing is intimately related to the thickness of the thermal boundary 
layers at the walls, it is then understandable that thinner boundary layers 
should yield a lower optimum spacing. 

Table 4.18 also shows that, for fixed Gr^, the maximum Nuav and the op- 
timum Ar both decrease with an increasing buoyancy parameter. The cause 
of the decrease in the Nuav is clear: when Gth is kept fixed, an increase 
in the buoyancy parameter implies decreased Ren, which in turn implies 
greater -viscosity and, as Pr is fixed, a greater a. This in turn implies thick- 
ened boundary layers and hence lower temperature gradients and Nuav The 
thicker boundary layer also imply that the optimum spacing should increase 
(and the optimum Ar should decrease). 
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Figure 4.22: Variation of average Nusselt number with aspect ratio and buoy- 
ancy parameter, Gth/Rch'^, at Grg = 10“^, for the case of mixed convection. 
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Multiple Channel Case 

As was elucidated in Section 2.7, the optimum spacing for the multiple chan- 
nel case is the one that maximises the product NuavAr. However, in the 
mixed convection NuavAr reaches a maximum as Aj. tends to infinity in the 
limit of small plate-to-plate spacing. 

The reason for this is evident when we consider that Uo remains fixed for 
all inter-plate spacing. So the flow rate per unit width of the channel remains 
fixed. Yet when the spacing becomes small, the outflow temperature will tend 
to reach its maximum which is 0 = 1, i.e., (T = Tyj). So if we quantify the 
problem, the average heat transfer coefficient is. 



where, qyj is the local wall heat flux. 

For the symmetrically heated walls case, the overall rate of heat transfer 
rate to the channel, Qch is given by 


Qch — 2 / Qyjdx 

Jo 

Therefore, 


j, Qch 

^“""2(r,„-Tx>)i? 

As the plate-to-plate spacing, L tends to zero, the fluid within the chan- 
nel gets heated to almost wall temperature, Tyj. Therefore, 


rISn Q<^ Cp {Tyj Tfxi) 


where the mass flow rate, m, is given by 
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m = pU(jL 


Therefore, when the spacing L -4 0 

U _P LcpUo 
~ 2H 

The average Nusselt number, Nuav is then given by 

hajjH PtRbjj 

= — = -2Ar 

where k is the thermal conductivity of the fluid (air). 
Thus, 


lim Nuav = 

Ar-^OO 


PrReu 

2Ar 


(4.3) 


So in the limit of large Ar, NuavAr tends to become a known constant. 
This can serve as a check for the accuracy of our simulation in the mixed 
convection case. Table 4.19 shows the limiting values of NuavAr obtained in 
the computation compared with the corresponding theoretical values. The 
comparison is quite good, indicating the accuracy of the computation. 


Figures 4.24, 4.25, and 4.26 show the product NuavAr versus Ar for 
the various buoyancy parameter and for different Grashof numbers {Gth = 
10®, 10^, and 10®). The Table 4.20 shows the aspect ratio {Ar)min.opt beyond 
which the value of NuavAr falls within 3% of the maximum value. Thus, any 
aspect ratio which is equal to or greater than the given values of (Ar) min.opt 
in Table 4.20 would be close to optimum, from the point of view of heat 
transfer. However, an increase in aspect ratio would cause a greater penalty, 
due to the higher pressure drop in a narrower channel. So it would be better 
to keep the aspect ratio very close to {Ar)min.opt a practical situation. 
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Table 4.19: Comparison of the theoretical and numerical values of the max- 
imum product of average Nusselt number and aspect ratio for different 
Grashof numbers and buoyancy parameters 


Gvh 

Gr. 

(^N 

Theoretical 

Numerical 


0.1 

35.00 

34.27 

10^ 

0.5 

15.65 

15.62 


1.0 

11.07 

11.19 


1.5 

9.04 

9.55 


0.1 

110.68 

108.37 

o 

1 — 1 

0.5 

49.50 

48.05 


1.0 

35.00 

33.84 


1.5 

28.58 

27.58 


0.1 

350.00 

343.53 

10^ 

0.5 

156.52 

153.69 


1.0 

110.68 

108.39 


1.5 

90.37 

88.35 


Table 4.20: Optimum aspect ratio for different Grashof numbers and buoy- 
ancy parameters: multiple channel case 


Grn 

min. opt 

II 

pi 

f-^i 

i 

P^ = 0.5 


= 1-5 


6.67 

5.00 

4.00 

3.33 

10^ 

10.00 

6.67 

5.56 

5.00 

10® 

25.00 

12.50 

10.00 

8.33 
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4.2.3 Pressure Drop 

In the case of mixed convection there is a pressure drop expected between 
the inlet and the outlet of the channel. The numerical computation allow 
the pressure drop to be calculated by a balance between the viscous, inertial 
and buoyancy forces using 


CD 

II 

<1 


dY - 


\/GrH 1 

dU 


Rejj [1 

0 W 


U=o dY + 0 dXdY 

dU 1 


+ 


where, APa„ is the nondimensionalised pressure (= . The value 

X = 1 + Cf appearing in the integral limits denotes the entrance to the 
computational domain which extends 10% of the channel height below the 
plates (i.e., ei =0.1). The pressure drop, obtained for various values of the 
buoyancy parameter and Gth = 10^, 10^, and 10® versus Ar is presented in 
Figure 4.27, 4.28, and 4.29. In all cases, it can be seen that the pressure 
drop is substantially less than what could be expected in the fully developed 
forced convection limit (which would be APav = Thus, buoyancy plays 
a very important role in reducing the pressure drop as compared to the forced 
convection case of equivalent Ren- 


The trends, for a given Gth and buoyancy parameter, in Figures 4.27, 
4.28, and 4.29 are easily understandable. The pressure drop increases with 
the aspect ratio as a narrower channel would, obviously, cause an increased 
pressure drop. For a given Gth and fixed Ar, the APa„ typically rises with 
increase in buoyancy parameter. This is because higher buoyancy parameter 
for fixed Gtjj corresponds to lower Ren and hence greater frictional loss. 

However, for some low aspect ratios, APa„ actually decreases when the 
buoyancy parameter is increased firom 0.5 to 1.0 and also from 1.0 to 1.5 in 
ail the figures. This may seem to be anomalous. However, this too can be 
explained: at low aspect ratios the channel spacing is wide and the thermal 
boundary layers of the opposite walls do not meet within the channel. The 
thermal boundary layers for higher buoyancy parameters would be thicker 
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Figure 4.24; Variation of the product of average Nusselt number and aspect 
ratio with aspect ratio and buoyancy parameter, Gru/Ren^, at Grii = 10®, 
for the case of mixed convection. 
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Figure 4.25: Variation of the product of average Nusselt number and aspect 
ratio with aspect ratio and buoyancy parameter, Gru/RcH^, at Gth = 
for the case of mixed convection. 
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0 5 10 15 20 25 30 35 

A, 


Figure 4.26: Variation of the product of average Nusselt number and aspect 
ratio with aspect ratio and buoyancy parameter, GrfffReH^, at Gth = 10®) 
for the case of mixed convection. 
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and, therefore, heat more of the channel fluid. The enhanced buoyancy ef- 
fect due to this seems to more than compensate for the increased viscosity 
(i.e. lower Reynolds number) at higher buoyancy parameter values. This in- 
creased buoyancy leads to a decrease in the pressure drop required to impose 
the flow. However, this effect does not carry over for higher aspect ratios 
when the thermal boundary layers meet within the channel and the major 
portion of the channel fluid is anyway heated for all buoyancy parameter 
values. 


4.2.4 Flow Rates in Natural and Mixed Convection 

The flow rate, for the natural convection case is compared with the mixed 
convection cases for various buoyancy parameters and for all Gth = 10^, 10^, 
and 10® in Figures 4.30, 4.31, and 4.32. In all cases, the is nondiraension- 
alised as in the natural convection case for the sole purpose of comparison 
(i.e., A4f = , where v is the volume flow rate per unit channel depth). 

It can be seen that the mixed convection flow rates are higher than in the 
natural convection cases. This is expected, as the pressure drop is positive 
(i.e., the outlet pressure is lower than the inlet pressure) in the mixed con- 
vection cases, as shown in Figures 4.27, 4.28, and 4.29, compared to the zero 
pressure drop in the natural convection cases. The pressure thus partially 
forces the flow in mixed convection, yielding a higher volume flow rate. 
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Figure 4.27: Variation of average pressure drop with aspect ratio and buoy- 
ancy parameter, GrulRen'^, at Gm = 10^, for the case of mixed convection. 
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0 2 4 6 8 10 12 14 16 18 20 

Ar 

Figure 4.28: Variation of average pr^sure drop with aspect ratio and buoy- 
ancy parameter, Grn / Ren^ , at Gr^ = 10^, for the case of mixed convection. 







Figure 4.30: Variation of imposed volume flow rate with aspect ratio and 
buoyancy parameter, Grji/ReH^, at Gr^ = 10^ and their comparison with 
the induced volume flow rate in the natural convection case. 
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Figure 4.31: Variation of imposed volume flow rate with aspect ratio and 
buoyancy parameter, Gru/Reff"^, at Crg = 10^ and their comparison with 
the induced volume flow rate in the natural convection case. 
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Figure 4.32: Variation of imposed volume flow rate with aspect ratio and 
buoyancy parameter, CrnlReH^, at Grn = 10® and their comparison with 
the induced volume flow rate in the natural convection case. 
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Chapter 5 

Conclusions and Scope for 
Future Work 


5.1 Conclusions 

The determination of optimum spacing between heated vertical parallel plates 
mounted equidistantly in an array of specified volume is a fundamental ther- 
mal design problem associated with the efficient cooling o nne ea ex 
change surfaces, nuclear reactors and electronic packages. Despite these im- 
portant applications, this problem had not been addressed for mixed convec- 
tion cooling. Although a number of studies of this problem with natural con- 
vection cooling have been reported in the literature, the volume^nstrained 
optimisation of inter-plate spacing has not received attention. These unex- 
amined aspects associated with optimisation of mter-plate spacing consti u 
the central motivation for the present investigation. 

To simulate the air cooling of electronic packages, an array of vertical, 
parallel, equidistant electronic circuit boards with el^tronic componen 
mounted on both sides of each board has been rnodeU^ as 
metrically heated, vertical, equal-spaced smooth parallel plates ^arntmu^ 
uniformly at constant temperature. There can exist ^o different optiimm 
inter-plate spacings such that the overall heat tra^ er is maximum. _ ^ 
corresponds to plates arranged in a constrained volume whereas ^ 
is associated with plates in an unconstrained volume. The latter optimu 
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spacing is essentially applicable to a single channel formed by a pair of sym- 
metrically heated vertical parallel plates. Consideration is given in this thesis 
to both constrained volume as well as unconstrained volume optimisation of 
inter-plate spacing. 

The optimum inter-plate spacing both for natural as well as mixed con- 
vection cooling have been determined by numerically obtaining the flow and 
thermal fields between two adjacent plates. The numerical results have been 
presented and discussed for a range of Grashof numbers (for natural convec- 
tion cooling), and combinations of Grashof number and buoyancy parameter 
(= GthIR^h) values (for mixed convection cooling), for a fluid Prandtl num- 
ber of 0.70. On the basis of the results, the following qualitative conclusions 
can drawn: 


1. The optimum spacing is intimately related to the boundary layer thick- 
ness at the two plates. The optimum spacing is typically reached at a 
configuration where the boundary layers meet well within the channel. 

2. The volume constrained optimum spacing is lower that the volume 
unconstrained optimum, in all cases. 

3. In natural convection cooling, the volume constrained optimum is stronger 
(i.e., more clearly evident), while in mixed convection it is the volume 
unconstrained optimum that is so. 

4. In mixed convection, the volume constrained spacing attains an asymp- 
totic optimum value as the spacing tends to zero. Beyond a limit, de- 
creases in inter-plate spacing results in higher pressure-drop without 
significant increase in overall heat transfer. 

5. The weak optimum at finite spacing (rather than an asymptotic op- 
timum at infinity) seen in natural convection is possibly due to the 
‘chimney effect’. 

6. For natural convection cooling, the optimum spacing decreases with 
increasing Grashof number. This is true for both volume constrained 
as well as volume unconstrained optimisation. The maximum Nusselt 
numbers attained at the optima are strongly and positively correlated 
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with the Grashof number. 


7. For mixed convection cooling, both optimum spacings decrease with 
increasing Grashof number, for a fixed buoyancy parameter value. The 
maximum Nusselt numbers attained are strongly and positively corre- 
lated with the Grashof number. 


8. For mixed convection cooling, both optimum spacings increase with in- 
creasing buoyancy parameter values, for a fibced Grashof number. The 
maximum Nusselt numbers attained are strongly and inversely corre- 
lated with the buoyancy parameter value (and thus strongly and posi- 
tively correlated with the Re 3 molds number). 


The quantitative results for specific parameter values are in Chapter 4. 


5.2 Scope of Future Work 

The aspects to be considered in the study of heated parallel plate channels 
are many, and greatly depend on the physical model involved. In the present 
investigation, the mathematical model has been obtained on the basis of 
certain approximations and idealisations presented in Section 2.3. Although, 
the present study brings out many important features related to thermal 
optimisation of spacing between heated parallel plates, the deviation of the 
present model from realistic applications in electronic cooling needs to be 
examined in future work. In the light of this, future research may be along 
the following lines: 

1. The problem of optimum spacing between electronic circuit boards with 
mounted components may be further studied by modelling them as 
heated plates with either protruding tyi>e or flush-mounted discrete 
heat sources. This problem can be studied for both natural as well as 
mixed convection cooling, with either constant heat flux or constant 
temperature wall conditions. 



142 chapter 5 . CONCLUSIONS AND SCOPE FOR FUTURE WORK 

hartnTc electronic packages are composed of circuit boards 
‘raaic nacTr “’y o»e side of each board. Such ela 

P'ote <ianne!f ““elelled as a series of asymmetrically heated 

form wa^rnmoe ““-T ° um 

‘P be lat a td ,h wrth uniform heat flux, and the other wall 
“efwal and Problem of optimum spacing can be studied for 

copsidtw d r “P be e>=‘endrf t 

ing discrete heat sources, as suggested above. ^ 

3 In 

‘be beeted parallel plate channels may be mod- 
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channels. ‘b* b®^' transfer interactions with the neighbouring 
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investigated for such situations. 
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